
# # # # # # # # # # # # # # # # # # 
#
# Replication Script
# 
# Part of replication of:
# Traditional Authorities, Norm Collisions, and Communal Conflict 
#
# Journal of Conflict Resolution
#
# By Clara Neupert-Wentz, 2024
#
# # # # # # # # # # # # # # # # # # 




###### Set-up ----

# Libraries
library(lme4)
library(lmerTest)
library(stargazer)
library(mediation)
library(interplot)
library(sensemakr)
library(tibble) 
library(insight)
library(texreg)
library(ggpubr) 
library(patchwork) 
library(fixest)
library(tidyverse)
library(haven)

# Cleanup
rm(list = ls())

# Globals
set.seed(9999)

setwd("/Neupert-Wentz Replication/")
fig.path <- "Neupert-Wentz Replication/figures"
tab.path <- "Neupert-Wentz Replication/tables"


#Load replication data
load("data/nc_rep.RData")


# Main analysis: Subset the data to groups that have some form of traditional political organization
nc_tpi <- nc_rep %>% 
  filter(any_tpi==1)




####### Descriptives ######


## Descriptives end of section 4.3 
nc_tpi %>% 
  filter(year==2015)  %>% 
  summarise(mean = mean(policing_subgroups, na.rm=T),
            median = median(policing_subgroups, na.rm=T),
            sd = sd(policing_subgroups, na.rm=T)) 

nc_tpi_15 <- nc_tpi %>% 
  filter(year==2015)  
table(nc_tpi_15$policing_subgroups)




#### Table A1: Summary Statistics ----


fileConn <- file()
writeLines(
  stargazer(as.data.frame(nc_tpi[c("n_events_acled_log", "mean_neigh_acled","n_events_ucdp_log", "mean_neigh_ucdp",
                                   "n_events_issue_territory_log", "mean_neigh_issue_territory",
                                   "n_events_issue_authority_log", "mean_neigh_issue_authority",
                                   "n_event_nongeo_log",
                                   "policing_capacity", "disp_capacity", "sec_capacity", "policing_subgroups",
                                   "policing_conflict",
                                   "two_class01_mcv", "pc1_mcv",
                                   "REGA_n_mcv",
                                   "status_egip", "capdist_wmean_ex_ln",
                                   "apop_ln", "sqkm_ln", "nightlight_corr_ex_ln", 
                                   "vdem_egaldem", "colpast", "max", "incidence_flag",
                                   "sub_ct",
                                   "year", "cow")]),
            title="Descriptive statistics", digits=2,
            covariate.labels=c("Communal Conflict (ACLED)","Spatial Lag (ACLED)", "Communal Conflict (UCDP)","Spatial Lag (UCDP)",
                               "Territory Conflict (UCDP)", "Territory Spatial Lag (UCDP)", "Authority Conflict (UCDP)", "Authority Spatial Lag (UCDP)", 
                               "Non-Geo Communal Conflict (UCDP)",
                               "Supergroup Policing Capacity", "Dispute Resolution Capacity Only", "Security Capacity Only", 
                               "No. of Subgroups w/ Policing Institutions",
                               "Jurisdictional Conflict",
                               "Constitutional Regulation (LCA)","Constitutional Regulation (PCA)","Constitutional Regulation (Min-Max)",
                               "Status Included (EPR)", "Capital Distance (log)",
                               "Population Size (log)", "Settlement Area (log)", "Night Light Density (log)", 
                               "Democracy (VDem)", "Colonial Past", "TG Population Share within State", "Rebellion Incidence",
                               "Subgroup Count (AMAR)",
                               "Year", "Country"),
            omit.summary.stat = c("p25", "p75"), 
            summary = T), 
  fileConn)
close(fileConn)



 ##### Analysis Set-up  ------

# Variables
outcomes_lags <- list(c("n_events_acled_log", "mean_neigh_acled"), #1
                      c("n_events_ucdp_log", "mean_neigh_ucdp"), #2
                      c("n_events_issue_territory_log", "mean_neigh_issue_territory"), #3
                      c("n_events_issue_authority_log", "mean_neigh_issue_authority"), #4
                      c("n_events_acled_log", "mean_neigh_acled", "policing_conflict"), #5
                      c("n_events_ucdp_log", "mean_neigh_ucdp", "policing_conflict"), #6
                      c("n_events_acled_log", "mean_neigh_acled", "sub_ct"), #7
                      c("n_events_ucdp_log", "mean_neigh_ucdp", "sub_ct"), #8
                      c("n_event_nongeo_log"), #9
                      c("n_events_acled_log", "mean_neigh_acled", "any_tpi"), #10
                      c("n_events_ucdp_log", "mean_neigh_ucdp", "any_tpi"), #11
                      c("policing_conflict") #12
                      ) 

treat <- c("policing_capacity", "policing_subgroups")
treat.1 <- "policing_capacity"
treat.2 <- "policing_subgroups"
treat.inter <- "policing_capacity*policing_subgroups"

cross.inter.list <- list(
  c("factor(two_class01_mcv)*policing_capacity", "policing_subgroups"), #1
  c("factor(two_class01_mcv)*policing_subgroups", "policing_capacity"), #2
  c("factor(two_class01_mcv)/policing_capacity", "policing_subgroups"), #3
  c("factor(two_class01_mcv)/policing_subgroups", "policing_capacity"), #4
  c("pc1_mcv*policing_capacity", "policing_subgroups"), #5
  c("pc1_mcv*policing_subgroups", "policing_capacity"), #6
  c("REGA_n_mcv*policing_capacity", "policing_subgroups"), #7
  c("REGA_n_mcv*policing_subgroups", "policing_capacity") #8
)


baseline.twfe <- "year + cow"
controls.twfe <- c("status_egip", "capdist_wmean_ex_ln",
                   "apop_ln", "sqkm_ln", "nightlight_corr_ex_ln", 
                   "vdem_egaldem", "incidence_flag")


baseline.me <- c("(1|year)", "(1|cow)")
controls.me <- c("status_egip", "capdist_wmean_ex_ln",
                 "apop_ln", "sqkm_ln", "nightlight_corr_ex_ln", 
                 "vdem_egaldem", "colpast", "max", "incidence_flag")


# Create formulas

# Main Effects
construct_formula <- function(outcome, treat, baseline, controls = NULL) {
  response = outcome[1]
  if(length(outcome) == 1){
    predictors = paste(treat, collapse = " + ")
  }
  if(length(outcome) == 2){
    predictors = paste(c(treat, outcome[2]), collapse = " + ")
  }
  if(length(outcome) == 3){
    predictors = paste(c(treat, outcome[2], outcome[3]), collapse = " + ")
  }
  if (!is.null(controls)) {
    predictors <- paste(predictors, paste(controls, collapse = " + "), sep = " + ")
  }
  formula = paste0(response, " ~ ", predictors, "|", baseline)
  as.formula(formula)
}

formulae <- list()

# Use lists to generate the formulas
formulae$twfe$twfe.1 <- lapply(outcomes_lags, construct_formula, treat = treat.1, baseline = baseline.twfe)
formulae$twfe$twfe.2 <- lapply(outcomes_lags, construct_formula, treat = treat.2, baseline = baseline.twfe)
formulae$twfe$twfe.3 <- lapply(outcomes_lags, construct_formula, treat = treat, baseline = baseline.twfe)
formulae$twfe$twfe.4 <- lapply(outcomes_lags, construct_formula, treat = treat, baseline = baseline.twfe, controls = controls.twfe)
formulae$twfe$twfe.5 <- lapply(outcomes_lags, construct_formula, treat = treat.inter, baseline = baseline.twfe, controls = controls.twfe)


# Cross-Level Interactions

construct_formula_inter <- function(inter, x) {
  if(length(x) == 1){
    return(reformulate(c(inter, baseline.me), response = x[1]))
  }
  if(length(x) == 2){
    return(reformulate(c(inter, x[2], baseline.me), response = x[1]))
  }
  if(length(x) == 3){
    return(reformulate(c(inter, x[2], x[3], baseline.me), response = x[1]))
  }
}

construct_formula_inter_controls <- function(inter, x) {
  if(length(x) == 1){
    return(reformulate(c(inter, baseline.me, controls.me), response = x[1]))
  }
  if(length(x) == 2){
    return(reformulate(c(inter, x[2], baseline.me, controls.me), response = x[1]))
  }
  if(length(x) == 3){
    return(reformulate(c(inter, x[2], x[3], baseline.me, controls.me), response = x[1]))
  }
}

formulae.me.inter.1 <- sapply(cross.inter.list, function(inter) {
  sapply(outcomes_lags, construct_formula_inter, inter = inter)
})


formulae.me.inter.2 <- sapply(cross.inter.list, function(inter) {
  sapply(outcomes_lags, construct_formula_inter_controls, inter = inter)
})




#### Main Estimation and Results ----


## > Figure 2 ----

# ACLED
model.ls <- list(
  lm.acled.1 <- feols(formulae$twfe$twfe.1[[1]],
                        data = nc_tpi, vcov = "hetero"),
  lm.acled.2 <- feols(formulae$twfe$twfe.2[[1]],
                        data = nc_tpi, vcov = "hetero"),
  lm.acled.3 <- feols(formulae$twfe$twfe.3[[1]],
                        data = nc_tpi, vcov = "hetero"),
  lm.acled.4 <- feols(formulae$twfe$twfe.4[[1]],
                        data = nc_tpi, vcov = "hetero")
) 


#Extract estimates
estimates.main_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Main", "Main","Main","Main","Main","Main"),
  estimate = c(coef(lm.acled.2)["policing_subgroups"],
               coef(lm.acled.3)["policing_subgroups"],
               coef(lm.acled.4)["policing_subgroups"],
               coef(lm.acled.1)["policing_capacity"],
               coef(lm.acled.3)["policing_capacity"],
               coef(lm.acled.4)["policing_capacity"]),
  ci.lower = c(confint(lm.acled.2)["policing_subgroups", 1],
               confint(lm.acled.3)["policing_subgroups", 1],
               confint(lm.acled.4)["policing_subgroups", 1],
               confint(lm.acled.1)["policing_capacity", 1],
               confint(lm.acled.3)["policing_capacity", 1],
               confint(lm.acled.4)["policing_capacity", 1]),
  ci.upper = c(confint(lm.acled.2)["policing_subgroups", 2],
               confint(lm.acled.3)["policing_subgroups", 2],
               confint(lm.acled.4)["policing_subgroups", 2],
               confint(lm.acled.1)["policing_capacity", 2],
               confint(lm.acled.3)["policing_capacity", 2],
               confint(lm.acled.4)["policing_capacity", 2])   )   



# UCDP
model.ls <- list(
  lm.ucdp.1 <- feols(formulae$twfe$twfe.1[[2]],
                data = nc_tpi, vcov = "hetero"),
  lm.ucdp.2 <- feols(formulae$twfe$twfe.2[[2]],
                data = nc_tpi, vcov = "hetero"),
  lm.ucdp.3 <- feols(formulae$twfe$twfe.3[[2]],
                data = nc_tpi, vcov = "hetero"),
  lm.ucdp.4 <- feols(formulae$twfe$twfe.4[[2]],
                data = nc_tpi, vcov = "hetero")
)


estimates.main_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Main", "Main","Main","Main","Main","Main"),
  estimate = c(coef(lm.ucdp.2)["policing_subgroups"],
               coef(lm.ucdp.3)["policing_subgroups"],
               coef(lm.ucdp.4)["policing_subgroups"],
               coef(lm.ucdp.1)["policing_capacity"],
               coef(lm.ucdp.3)["policing_capacity"],
               coef(lm.ucdp.4)["policing_capacity"]),
  ci.lower = c(confint(lm.ucdp.2)["policing_subgroups", 1],
               confint(lm.ucdp.3)["policing_subgroups", 1],
               confint(lm.ucdp.4)["policing_subgroups", 1],
               confint(lm.ucdp.1)["policing_capacity", 1],
               confint(lm.ucdp.3)["policing_capacity", 1],
               confint(lm.ucdp.4)["policing_capacity", 1]),
  ci.upper = c(confint(lm.ucdp.2)["policing_subgroups", 2],
               confint(lm.ucdp.3)["policing_subgroups", 2],
               confint(lm.ucdp.4)["policing_subgroups", 2],
               confint(lm.ucdp.1)["policing_capacity", 2],
               confint(lm.ucdp.3)["policing_capacity", 2],
               confint(lm.ucdp.4)["policing_capacity", 2])   )  



estimates <- bind_rows(
  list(estimates.main_ucdp,
       estimates.main_acled)
)

estimates$model_factor <- factor(estimates$model, 
                                 levels = c("Incl. controls", "Both treatments", "One treatment"))

save(estimates, file="estimates/estimates.RData")

estimates %>%
  ggplot(aes(x=estimate, y=variable, shape=model_factor)) +
  geom_point(position = position_dodge(width=.3))+ 
  geom_errorbarh(aes(xmin=ci.lower, xmax=ci.upper), position=position_dodge(width=.3),
                 height=0) + 
  labs(x="", y="", shape="Specification:") + 
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_grid(. ~ test, scales = "free_x") +
  theme(panel.spacing = unit(2, "lines")) +
  scale_shape_discrete(breaks=c("One treatment", "Both treatments","Incl. controls")) 


ggsave("coef_plot_main.pdf", 
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =10,
       units = "cm",
       dpi=700)



## > Table A3 ----

etable(list(lm.acled.1, lm.acled.2, lm.acled.3, lm.acled.4, lm.ucdp.1, lm.ucdp.2, lm.ucdp.3, lm.ucdp.4),
       title = "Traditional Policing Institutions and Communal Conflict",
       dict = c(policing_capacity="Supergroup Policing", policing_subgroups="No. of Subgroups w/ Policing Institutions", 
                status_egip="\\rule{0pt}{3ex}Status Included (EPR)", 
                capdist_wmean_ex_ln="Capital Distance (log)", apop_ln="Population Size (log)",
                sqkm_ln="Settlement Area (log)", nightlight_corr_ex_ln="Night Light Density (log)",
                vdem_egaldem = "Democracy (VDem)", incidence_flag="Rebellion Incidence",
                 mean_neigh_acled="Spatial Lag (ACLED)", mean_neigh_ucdp="Spatial Lag (UCDP)"),
       order = c("%policing_capacity", "%policing_subgroups", "%status_egip", "%capdist_wmean_ex_ln", 
                 "%apop_ln", "%sqkm_ln", "%nightlight_corr_ex_ln", "%vdem_egaldem",
                 "%incidence_flag", "%mean_neigh_acled", "%mean_neigh_ucdp"),
       digits = "r3",
       notes = "Linear models.",
       #file = file.path("tables/main.tex"),
       replace = TRUE,
       label = "tab:main",
       style.tex = style.tex(adjustbox = TRUE))



##  > Figure 3 ----

# ACLED

#Interaction
model.ls <- list(
  inter.lmer.acled.1 <- lmer(formulae.me.inter.1[1, 1][[1]], 
                           data = nc_tpi),
  inter.lmer.acled.2 <- lmer(formulae.me.inter.2[1, 1][[1]], 
                           data = nc_tpi),
  inter.lmer.acled.3 <- lmer(formulae.me.inter.1[1, 2][[1]], 
                           data = nc_tpi),
  inter.lmer.acled.4 <- lmer(formulae.me.inter.2[1, 2][[1]], 
                           data = nc_tpi)
)

#Regulated vs. Unregulated
model.ls <- list(
  inter.lmer.acled.diff.1 <- lmer(formulae.me.inter.1[1, 3][[1]], 
                           data = nc_tpi),
  inter.lmer.acled.diff.2 <- lmer(formulae.me.inter.2[1, 3][[1]], 
                           data = nc_tpi),
  inter.lmer.acled.diff.3 <- lmer(formulae.me.inter.1[1, 4][[1]], 
                           data = nc_tpi),
  inter.lmer.acled.diff.4 <- lmer(formulae.me.inter.2[1, 4][[1]], 
                           data = nc_tpi)
)


# UCDP

#Interaction
model.ls <- list(
  inter.lmer.ucdp.1 <- lmer(formulae.me.inter.1[2, 1][[1]], 
                          data = nc_tpi),
  inter.lmer.ucdp.2 <- lmer(formulae.me.inter.2[2, 1][[1]], 
                          data = nc_tpi),
  inter.lmer.ucdp.3 <- lmer(formulae.me.inter.1[2, 2][[1]], 
                          data = nc_tpi),
  inter.lmer.ucdp.4 <- lmer(formulae.me.inter.2[2, 2][[1]], 
                          data = nc_tpi)
)

#Regulated vs. Unregulated
model.ls <- list(
  inter.lmer.ucdp.diff.1 <- lmer(formulae.me.inter.1[2, 3][[1]], 
                          data = nc_tpi),
  inter.lmer.ucdp.diff.2 <- lmer(formulae.me.inter.2[2, 3][[1]], 
                          data = nc_tpi),
  inter.lmer.ucdp.diff.3 <- lmer(formulae.me.inter.1[2, 4][[1]], 
                          data = nc_tpi),
  inter.lmer.ucdp.diff.4 <- lmer(formulae.me.inter.2[2, 4][[1]], 
                          data = nc_tpi)
)


#Extract estimates
estimates.inter <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("Main", "Main","Main","Main","Main","Main","Main","Main"),
  estimate = c(coef(inter.lmer.acled.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.acled.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.acled.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(inter.lmer.acled.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(inter.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.ucdp.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(inter.lmer.ucdp.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )


estimates.inter$model_factor <- factor(estimates.inter$model, 
                                       levels = c("Incl. controls", "Excl. controls"))

save(estimates.inter, file="estimates/estimates.inter.RData" )

#Extract estimates regulated vs. unregulated
estimates.inter.diff <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("Main", "Main","Main","Main","Main","Main","Main","Main"),
  estimate = c(coef(inter.lmer.acled.diff.1)$cow[1,"factor(two_class01_mcv)0:policing_capacity"],
               coef(inter.lmer.acled.diff.2)$cow[1,"factor(two_class01_mcv)0:policing_capacity"],
               coef(inter.lmer.acled.diff.3)$cow[1,"factor(two_class01_mcv)0:policing_subgroups"],
               coef(inter.lmer.acled.diff.4)$cow[1,"factor(two_class01_mcv)0:policing_subgroups"],
               coef(inter.lmer.ucdp.diff.1)$cow[1,"factor(two_class01_mcv)0:policing_capacity"],
               coef(inter.lmer.ucdp.diff.2)$cow[1,"factor(two_class01_mcv)0:policing_capacity"],
               coef(inter.lmer.ucdp.diff.3)$cow[1,"factor(two_class01_mcv)0:policing_subgroups"],
               coef(inter.lmer.ucdp.diff.4)$cow[1,"factor(two_class01_mcv)0:policing_subgroups"],
               coef(inter.lmer.acled.diff.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.acled.diff.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.acled.diff.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(inter.lmer.acled.diff.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(inter.lmer.ucdp.diff.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.ucdp.diff.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.ucdp.diff.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(inter.lmer.ucdp.diff.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(inter.lmer.acled.diff.1)["factor(two_class01_mcv)0:policing_capacity", 1],
               confint(inter.lmer.acled.diff.2)["factor(two_class01_mcv)0:policing_capacity", 1],
               confint(inter.lmer.acled.diff.3)["factor(two_class01_mcv)0:policing_subgroups", 1],
               confint(inter.lmer.acled.diff.4)["factor(two_class01_mcv)0:policing_subgroups", 1],
               confint(inter.lmer.ucdp.diff.1)["factor(two_class01_mcv)0:policing_capacity", 1],
               confint(inter.lmer.ucdp.diff.2)["factor(two_class01_mcv)0:policing_capacity", 1],
               confint(inter.lmer.ucdp.diff.3)["factor(two_class01_mcv)0:policing_subgroups", 1],
               confint(inter.lmer.ucdp.diff.4)["factor(two_class01_mcv)0:policing_subgroups", 1],
               confint(inter.lmer.acled.diff.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.acled.diff.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.acled.diff.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(inter.lmer.acled.diff.4)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(inter.lmer.ucdp.diff.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.ucdp.diff.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.ucdp.diff.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(inter.lmer.ucdp.diff.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(inter.lmer.acled.diff.1)["factor(two_class01_mcv)0:policing_capacity", 2],
               confint(inter.lmer.acled.diff.2)["factor(two_class01_mcv)0:policing_capacity", 2],
               confint(inter.lmer.acled.diff.3)["factor(two_class01_mcv)0:policing_subgroups", 2],
               confint(inter.lmer.acled.diff.4)["factor(two_class01_mcv)0:policing_subgroups", 2],
               confint(inter.lmer.ucdp.diff.1)["factor(two_class01_mcv)0:policing_capacity", 2],
               confint(inter.lmer.ucdp.diff.2)["factor(two_class01_mcv)0:policing_capacity", 2],
               confint(inter.lmer.ucdp.diff.3)["factor(two_class01_mcv)0:policing_subgroups", 2],
               confint(inter.lmer.ucdp.diff.4)["factor(two_class01_mcv)0:policing_subgroups", 2],
               confint(inter.lmer.acled.diff.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.acled.diff.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.acled.diff.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(inter.lmer.acled.diff.4)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(inter.lmer.ucdp.diff.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.ucdp.diff.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.ucdp.diff.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(inter.lmer.ucdp.diff.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )

estimates.inter.diff$model_factor <- factor(estimates.inter.diff$model, 
                                       levels = c("Incl. controls", "Excl. controls"))


save(estimates.inter.diff, file="estimates/estimates.inter.diff.RData" ) 


# Combine data frames
estimates.combined <- rbind(
  transform(estimates.inter, source = "Interaction"),
  transform(estimates.inter.diff, source =  c("Marginal Effect, Unregulated", "Marginal Effect, Unregulated","Marginal Effect, Unregulated","Marginal Effect, Unregulated",
             "Marginal Effect, Unregulated", "Marginal Effect, Unregulated","Marginal Effect, Unregulated","Marginal Effect, Unregulated",
             "Marginal Effect, Regulated","Marginal Effect, Regulated","Marginal Effect, Regulated","Marginal Effect, Regulated",
             "Marginal Effect, Regulated","Marginal Effect, Regulated","Marginal Effect, Regulated","Marginal Effect, Regulated"))
)

# Specify the levels
estimates.combined$source <- factor(estimates.combined$source, 
                                    levels = c("Marginal Effect, Regulated","Marginal Effect, Unregulated",  "Interaction"))

save(estimates.inter.combined, file="estimates/estimates.inter.combined.RData" )

# Color palette
color_palette <- c("Interaction" = "black", "Marginal Effect, Unregulated" = "gray", "Marginal Effect, Regulated" = "slategray")

# Linetype palette 
linetype_palette <- c("Interaction" = "solid", "Marginal Effect, Unregulated" = "solid", "Marginal Effect, Regulated" = "solid")

# Positioning 
dodge <- position_dodge(width=0.2)

# Plot
estimates.combined %>%
  ggplot(aes(x = estimate, y = variable, color = source, shape = model_factor, linetype = source)) +
  geom_point(position = dodge) +
  geom_errorbarh(aes(xmin = ci.lower, xmax = ci.upper), position = dodge, height = 0) +
  labs(x = "", y = "", shape = "Specification:", color = "Effects:", linetype = "Effects:") +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_grid(. ~ test, scales = "free_x") +
  theme(panel.spacing = unit(2, "lines")) +
  scale_shape_discrete(breaks = c("Excl. controls", "Incl. controls")) +
  scale_color_manual(values = color_palette) +
  scale_linetype_manual(values = linetype_palette) +
  guides(shape = guide_legend(order = 2)) 


ggsave("coef_plot_inter.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 21,
       height =10,
       units = "cm",
       dpi=700)



##  > Figure A4 ----
texreg(list(inter.lmer.acled.1, inter.lmer.acled.2, inter.lmer.acled.3, inter.lmer.acled.4, inter.lmer.ucdp.1, inter.lmer.ucdp.2, inter.lmer.ucdp.3, inter.lmer.ucdp.4), include.ci = FALSE,
       #file = "tables/inter.tex",
       stars = c(0.001, 0.01, 0.05, 0.1),
       label = "tab:inter",
       custom.coef.names = c("Intercept", "\\rule{0pt}{3ex}Constitutional Regulation",
                             "\\rule{0pt}{3ex}Supergroup Policing", 
                             "No. of Subgroups w/ Policing Institutions",
                             "Spatial Lag (ACLED)", 
                             "\\rule{0pt}{3ex}Supergroup Policing $\\times$  \\hspace{3mm}Constitutional Regulation",
                             "Status Included (EPR)",
                             "Capital Distance (log)", "Population Size (log)", "Settlement Area (log)", 
                             "Night Light Density (log)", "Democracy (VDem)", 
                             "Colonial Past", "TG Population Share within State",
                             "Rebellion Incidence",
                             "\\rule{0pt}{3ex}No. of Subgroups $\\times$  \\hspace{3mm}Constitutional Regulation",
                             "Spatial Lag (UCDP)"),
       reorder.coef = c(6, 16, 3, 4, 2, 7, 8, 9, 10, 11, 12, 13, 14, 15, 5, 17, 1),
       center = T,
       caption = "Traditional Policing Institutions , Constitutional Regulation, and Communal Conflict",
       scriptsize=T,
       caption.above = T,
       siunitx = F, 
       custom.gof.rows = list("\\rule{0pt}{3ex}Country {\\scshape RE}" = c("\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark"),
                              "Year {\\scshape RE}" = c("\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark")),
       include.adjrs = FALSE, include.bic = FALSE, include.variance = T,
       custom.note = "\\textit{Notes:} Linear mixed-effect models with random intercepts for countries and years. Standard errors in parentheses. $^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001."
)


### Mechanism: Outcome Vertical Jurisdictional Conflict ----

## > Figure 4 ----

# Cross-section formulas main effect
baseline.twfe.confl <-    "cow"
formulae$twfe$twfe.confl.1 <- lapply(outcomes_lags, construct_formula, treat = treat.1, baseline = baseline.twfe.confl)
formulae$twfe$twfe.confl.2 <- lapply(outcomes_lags, construct_formula, treat = treat.2, baseline = baseline.twfe.confl)
formulae$twfe$twfe.confl.3 <- lapply(outcomes_lags, construct_formula, treat = treat, baseline = baseline.twfe.confl)
formulae$twfe$twfe.confl.4 <- lapply(outcomes_lags, construct_formula, treat = treat, baseline = baseline.twfe.confl, controls = controls.twfe)


# Cross-section interaction formulas
baseline.me.confl   <-    "(1 | cow)" 

construct_formula_inter_confl <- function(inter, x) {
  if(length(x) == 1){
    return(reformulate(c(inter, baseline.me.confl), response = x[1]))
  }
  if(length(x) == 2){
    return(reformulate(c(inter, x[2], baseline.me.confl), response = x[1]))
  }
  if(length(x) == 3){
    return(reformulate(c(inter, x[2], x[3], baseline.me.confl), response = x[1]))
  }
}

construct_formula_inter_controls_confl <- function(inter, x) {
  if(length(x) == 1){
    return(reformulate(c(inter, baseline.me.confl, controls.me), response = x[1]))
  }
  if(length(x) == 2){
    return(reformulate(c(inter, x[2], baseline.me.confl, controls.me), response = x[1]))
  }
  if(length(x) == 3){
    return(reformulate(c(inter, x[2], x[3], baseline.me.confl, controls.me), response = x[1]))
  }
}

formulae.me.inter.confl.1 <- sapply(cross.inter.list, function(inter) {
  sapply(outcomes_lags, construct_formula_inter_confl, inter = inter)
})

formulae.me.inter.confl.2 <- sapply(cross.inter.list, function(inter) {
  sapply(outcomes_lags, construct_formula_inter_controls_confl, inter = inter)
})



# Models
model.ls <- list(
  lm.confl.1 <- feols(formulae$twfe$twfe.confl.1[[12]],
                      data = nc_tpi[nc_tpi$year == 2015,]),
  lm.confl.2 <- feols(formulae$twfe$twfe.confl.2[[12]],
                      data = nc_tpi[nc_tpi$year == 2015,]),
  lm.confl.3 <- feols(formulae$twfe$twfe.confl.3[[12]],
                      data = nc_tpi[nc_tpi$year == 2015,]),
  lm.confl.4 <- feols(formulae$twfe$twfe.confl.4[[12]],
                      data = nc_tpi[nc_tpi$year == 2015,])
)


model.ls <- list(
  inter.lmer.confl.1 <- lmer(formulae.me.inter.confl.1[12, 1][[1]], 
                             data = nc_tpi[nc_tpi$year == 2015,]),
  inter.lmer.confl.2 <- lmer(formulae.me.inter.confl.2[12, 1][[1]], 
                             data = nc_tpi[nc_tpi$year == 2015,]),
  inter.lmer.confl.3 <- lmer(formulae.me.inter.confl.1[12, 2][[1]], 
                             data = nc_tpi[nc_tpi$year == 2015,]),
  inter.lmer.confl.4 <- lmer(formulae.me.inter.confl.2[12, 2][[1]], 
                             data = nc_tpi[nc_tpi$year == 2015,])
)


#Extract estimates

estimates.main_confl <- data.frame(
  test = c( "Jurisdictional Conflict", "Jurisdictional Conflict", "Jurisdictional Conflict",
            "Jurisdictional Conflict", "Jurisdictional Conflict", "Jurisdictional Conflict"),
  model = c("One treatment", "One treatment", "Both treatments",  "Both treatments",
            "Incl. controls", "Incl. controls"),
  variable = c( "Supergroup Policing",   "No. of Subgroups \n w/ Policing Institutions", 
                "Supergroup Policing",  "No. of Subgroups \n w/ Policing Institutions", 
                "Supergroup Policing", "No. of Subgroups \n w/ Policing Institutions"),
  estimate = c(coef(lm.confl.1)["policing_capacity"],
               coef(lm.confl.2)["policing_subgroups"],
               coef(lm.confl.3)["policing_capacity"],
               coef(lm.confl.3)["policing_subgroups"],
               coef(lm.confl.4)["policing_capacity"],
               coef(lm.confl.4)["policing_subgroups"]),
  ci.lower = c(confint(lm.confl.1)["policing_capacity", 1],
               confint(lm.confl.2)["policing_subgroups", 1],
               confint(lm.confl.3)["policing_capacity", 1],
               confint(lm.confl.3)["policing_subgroups", 1],
               confint(lm.confl.4)["policing_capacity", 1],
               confint(lm.confl.4)["policing_subgroups", 1]),
  ci.upper = c(confint(lm.confl.1)["policing_capacity", 2],
               confint(lm.confl.2)["policing_subgroups", 2],
               confint(lm.confl.3)["policing_capacity", 2],
               confint(lm.confl.3)["policing_subgroups", 2],
               confint(lm.confl.4)["policing_capacity", 2],
               confint(lm.confl.4)["policing_subgroups", 2])   )   


#Extract estimates interactions

estimates.inter.confl <- data.frame(
  test = c("Jurisdictional Conflict", "Jurisdictional Conflict", "Jurisdictional Conflict", 
           "Jurisdictional Conflict"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  estimate = c(coef(inter.lmer.confl.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.confl.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(inter.lmer.confl.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(inter.lmer.confl.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(inter.lmer.confl.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.confl.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(inter.lmer.confl.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(inter.lmer.confl.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(inter.lmer.confl.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.confl.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(inter.lmer.confl.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(inter.lmer.confl.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )

estimates.confl <- bind_rows(
  list(estimates.main_confl,
       estimates.inter.confl)
)

estimates.confl$model_factor <- factor(estimates.confl$model, 
                                       levels = c("Excl. controls", "Incl. controls", "Both treatments", "One treatment"))

estimates.confl$variable_factor <- factor(estimates.confl$variable, 
                                          levels = c("No. of Subgroups x \n Constitutional Regulation",
                                                     "Supergroup Policing x \n Constitutional Regulation", 
                                                     "No. of Subgroups \n w/ Policing Institutions", 
                                                     "Supergroup Policing"
                                          ))


save(estimates.confl, file="estimates/estimates.confl.RData" ) 

estimates.confl %>%
  ggplot(aes(x=estimate, y=variable_factor, shape=model_factor)) +
  geom_point(position = position_dodge(width=.3))+ 
  geom_errorbarh(aes(xmin=ci.lower, xmax=ci.upper), position=position_dodge(width=.3),
                 height=0) + 
  labs(x="", y="", shape="Specification:") + 
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_grid(. ~ test, scales = "free_x") +
  theme(panel.spacing = unit(2, "lines")) +
  scale_shape_discrete(breaks=c("One treatment", "Both treatments","Incl. controls", "Excl. controls")) +
  xlim(-0.8, 0.8)


ggsave("coef_plot_confl.pdf", 
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =10,
       units = "cm",
       dpi=700)


##  > Table A5 ----

texreg(list(lm.confl.1, lm.confl.2, lm.confl.3, lm.confl.4, inter.lmer.confl.1, inter.lmer.confl.2, inter.lmer.confl.3, inter.lmer.confl.4), include.ci = FALSE,
       #file = "tables/confl.tex",
       stars = c(0.001, 0.01, 0.05, 0.1),
       label = "tab:confl",
       custom.coef.names = c("Supergroup Policing", 
                             "No. of Subgroups w/ Policing Institutions",
                             "Status Included (EPR)",
                             "Capital Distance (log)", "Population Size (log)", "Settlement Area (log)", 
                             "Night Light Density (log)",
                             "Rebellion Incidence",
                             "Intercept", 
                             "Constitutional Regulation",
                             "\\rule{0pt}{3ex}Supergroup Policing $\\times$  \\hspace{3mm}Constitutional Regulation",
                             "Democracy (VDem)", 
                             "Colonial Past", "TG Population Share within State",
                             "\\rule{0pt}{3ex}No. of Subgroups $\\times$  \\hspace{3mm}Constitutional Regulation"),
       center = T,
       reorder.coef = c(1, 2, 11, 15, 10, 3, 4, 5, 6, 7, 8, 12, 13, 14, 9),
       caption = "Traditional Policing and Jurisdictional Conflict",
       scriptsize=T,
       caption.above = T,
       siunitx = F, 
       custom.gof.rows = list("Country {\\scshape FE}" = c("\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "$\\times$", "$\\times$", "$\\times$", "$\\times$"),
                              "Country {\\scshape RE}" = c("$\\times$", "$\\times$", "$\\times$", "$\\times$", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark")),
       include.adjrs = FALSE, include.bic = FALSE, include.variance = T,
       custom.note = "\\textit{Notes:} Linear models. Standard errors in parentheses $^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001."
)


## > Mediation I: Controlling for Jurisdictional Conflict (Figure A3) ----

# ACLED models
model.ls <- list(
  med.lm.acled.1 <- feols(formulae$twfe$twfe.1[[5]],
                          data = nc_tpi, vcov = "hetero"),
  med.lm.acled.2 <- feols(formulae$twfe$twfe.2[[5]],
                          data = nc_tpi, vcov = "hetero"),
  med.lm.acled.3 <- feols(formulae$twfe$twfe.3[[5]],
                          data = nc_tpi, vcov = "hetero"),
  med.lm.acled.4 <- feols(formulae$twfe$twfe.4[[5]],
                          data = nc_tpi, vcov = "hetero")
)

# Extract estimates
estimates.med_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Controlling for Jurisdictional Conflict", "Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict"),
  estimate = c(coef(med.lm.acled.2)["policing_subgroups"],
               coef(med.lm.acled.3)["policing_subgroups"],
               coef(med.lm.acled.4)["policing_subgroups"],
               coef(med.lm.acled.1)["policing_capacity"],
               coef(med.lm.acled.3)["policing_capacity"],
               coef(med.lm.acled.4)["policing_capacity"]),
  ci.lower = c(confint(med.lm.acled.2)["policing_subgroups", 1],
               confint(med.lm.acled.3)["policing_subgroups", 1],
               confint(med.lm.acled.4)["policing_subgroups", 1],
               confint(med.lm.acled.1)["policing_capacity", 1],
               confint(med.lm.acled.3)["policing_capacity", 1],
               confint(med.lm.acled.4)["policing_capacity", 1]),
  ci.upper = c(confint(med.lm.acled.2)["policing_subgroups", 2],
               confint(med.lm.acled.3)["policing_subgroups", 2],
               confint(med.lm.acled.4)["policing_subgroups", 2],
               confint(med.lm.acled.1)["policing_capacity", 2],
               confint(med.lm.acled.3)["policing_capacity", 2],
               confint(med.lm.acled.4)["policing_capacity", 2])   )   


# UCDP models
model.ls <- list(
  med.lm.ucdp.1 <- feols(formulae$twfe$twfe.1[[6]],
                         data = nc_tpi, vcov = "hetero"),
  med.lm.ucdp.2 <- feols(formulae$twfe$twfe.2[[6]],
                         data = nc_tpi, vcov = "hetero"),
  med.lm.ucdp.3 <- feols(formulae$twfe$twfe.3[[6]],
                         data = nc_tpi, vcov = "hetero"),
  med.lm.ucdp.4 <- feols(formulae$twfe$twfe.4[[6]],
                         data = nc_tpi, vcov = "hetero")
)

# Extract estimates
estimates.med_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Controlling for Jurisdictional Conflict", "Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict"),
  estimate = c(coef(med.lm.ucdp.2)["policing_subgroups"],
               coef(med.lm.ucdp.3)["policing_subgroups"],
               coef(med.lm.ucdp.4)["policing_subgroups"],
               coef(med.lm.ucdp.1)["policing_capacity"],
               coef(med.lm.ucdp.3)["policing_capacity"],
               coef(med.lm.ucdp.4)["policing_capacity"]),
  ci.lower = c(confint(med.lm.ucdp.2)["policing_subgroups", 1],
               confint(med.lm.ucdp.3)["policing_subgroups", 1],
               confint(med.lm.ucdp.4)["policing_subgroups", 1],
               confint(med.lm.ucdp.1)["policing_capacity", 1],
               confint(med.lm.ucdp.3)["policing_capacity", 1],
               confint(med.lm.ucdp.4)["policing_capacity", 1]),
  ci.upper = c(confint(med.lm.ucdp.2)["policing_subgroups", 2],
               confint(med.lm.ucdp.3)["policing_subgroups", 2],
               confint(med.lm.ucdp.4)["policing_subgroups", 2],
               confint(med.lm.ucdp.1)["policing_capacity", 2],
               confint(med.lm.ucdp.3)["policing_capacity", 2],
               confint(med.lm.ucdp.4)["policing_capacity", 2])   )  



### Interaction models

# ACLED
model.ls <- list(
  med.inter.lmer.acled.1 <- lmer(formulae.me.inter.1[5, 1][[1]], 
                                 data = nc_tpi),
  med.inter.lmer.acled.2 <- lmer(formulae.me.inter.2[5, 1][[1]], 
                                 data = nc_tpi),
  med.inter.lmer.acled.3 <- lmer(formulae.me.inter.1[5, 2][[1]], 
                                 data = nc_tpi),
  med.inter.lmer.acled.4 <- lmer(formulae.me.inter.2[5, 2][[1]], 
                                 data = nc_tpi)
)


# UCDP
model.ls <- list(
  med.inter.lmer.ucdp.1 <- lmer(formulae.me.inter.1[6, 1][[1]], 
                                data = nc_tpi),
  med.inter.lmer.ucdp.2 <- lmer(formulae.me.inter.2[6, 1][[1]], 
                                data = nc_tpi),
  med.inter.lmer.ucdp.3 <- lmer(formulae.me.inter.1[6, 2][[1]], 
                                data = nc_tpi),
  med.inter.lmer.ucdp.4 <- lmer(formulae.me.inter.2[6, 2][[1]], 
                                data = nc_tpi)
)



#Extract estimates
estimates.inter.med <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("Controlling for Jurisdictional Conflict", "Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict","Controlling for Jurisdictional Conflict"),
  estimate = c(coef(med.inter.lmer.acled.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(med.inter.lmer.acled.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(med.inter.lmer.acled.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(med.inter.lmer.acled.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(med.inter.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(med.inter.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(med.inter.lmer.ucdp.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(med.inter.lmer.ucdp.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(med.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(med.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(med.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(med.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(med.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(med.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(med.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(med.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(med.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(med.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(med.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(med.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(med.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(med.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(med.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(med.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )

estimates.med <- bind_rows(
  list(estimates.med_ucdp,
       estimates.med_acled,
       estimates.inter.med)
)


estimates.med$model_factor <- factor(estimates.med$model, 
                                     levels = c("Incl. controls", "Both treatments", "One treatment"))
# Presented in Figure A3
save(estimates.med, file="estimates/estimates.med.RData" )



####### > Mediation II: Figure 5 and Figure A10  ----

# ACLED
fit_data <- nc_tpi[!is.na(nc_tpi$n_events_acled_log),]
fit_data <- fit_data %>% select(n_events_acled_log, 
                                policing_capacity, cow , year ,
                                mean_neigh_acled ,  status_egip , capdist_wmean_ex_ln , apop_ln , sqkm_ln , 
                                nightlight_corr_ex_ln , vdem_egaldem , incidence_flag ,
                                policing_conflict)
fit_data <- fit_data[complete.cases(fit_data),]

# fixed-effect-demean data by hand, to solve all kinds of problems with lm_robust & mediate
fit_data <- demean(fit_data %>% select(-c(cow, year)), 
                   fit_data %>% select(cow, year))

Y_acled <- lm(n_events_acled_log ~  policing_capacity  + 
                mean_neigh_acled +  status_egip  + capdist_wmean_ex_ln + apop_ln + sqkm_ln + 
                nightlight_corr_ex_ln + vdem_egaldem + incidence_flag +
                policing_conflict,
                data= fit_data)

M_acled <- lm(policing_conflict ~  policing_capacity + 
                mean_neigh_acled +  status_egip  + capdist_wmean_ex_ln + apop_ln + sqkm_ln + 
                nightlight_corr_ex_ln + vdem_egaldem + incidence_flag,
              data= fit_data)

mediation.result_acled <- mediate(M_acled, Y_acled, sims=1000, treat="policing_capacity", 
                                  robustSE = T,
                                  mediator = "policing_conflict")

summary(mediation.result_acled)

# UCDP
fit_data <- nc_tpi[!is.na(nc_tpi$n_events_ucdp_log),]
fit_data <- fit_data %>% select(n_events_ucdp_log, 
                                policing_capacity, cow , year ,
                                mean_neigh_ucdp ,  status_egip  , capdist_wmean_ex_ln , apop_ln , sqkm_ln , 
                                nightlight_corr_ex_ln , vdem_egaldem , incidence_flag ,
                                policing_conflict)
fit_data <- fit_data[complete.cases(fit_data),]

# fixed-effect-demean data by hand, to solve all kinds of problems with lm_robust & mediate
fit_data <- demean(fit_data %>% select(-c(cow, year)), 
                   fit_data %>% select(cow, year))

Y_ucdp <- lm(n_events_ucdp_log ~  policing_capacity  + #factor(cow) + factor(year) +
               mean_neigh_ucdp +  status_egip  + capdist_wmean_ex_ln + apop_ln + sqkm_ln + 
               nightlight_corr_ex_ln + vdem_egaldem + incidence_flag +
               policing_conflict, data= fit_data)

M_ucdp <- lm(policing_conflict ~  policing_capacity   + #factor(cow) + factor(year)+
               mean_neigh_ucdp +  status_egip  + capdist_wmean_ex_ln + apop_ln + sqkm_ln + 
               nightlight_corr_ex_ln + vdem_egaldem + incidence_flag,
             data=fit_data)

mediation.result_ucdp <- mediate(M_ucdp, Y_ucdp, sims=1000, 
                                 robustSE = T,
                                 treat="policing_capacity", 
                                 mediator = "policing_conflict")

summary(mediation.result_ucdp)

pdf(file= "figures/mediation.pdf",  
    height=9, width=12)
par(mfrow = c(1, 2))
plot(mediation.result_acled,
     main = "ACLED")

plot(mediation.result_ucdp,
     main = "UCDP")
dev.off() 




##### Appendix G: Sensitivity Analysis -----


## > Figure A6: Sensitivity of Policing Capacity  ------

# Relevant models are specified above:
  # lm.acled.4
  # lm.ucdp.4

# function to compute bounds manually - for standard models
# (Direct effect models below need another function)
plot_sensitivity <- function(controlvar, treatvar,
                             model, labels, type, k = 3){
  
  for(control in controlvar){
    
    controlvar.index <- which(controlvar == control)
    
# Use the t statistic of control variable in the outcome regression
# to compute the partial R2 of control variable with the outcome.
    if(model[controlvar.index] == "lm.acled.4"){
      dof <- degrees_freedom(lm.acled.4, type = "resid")
      r2yxj.dx <- partial_r2(t_statistic = lm.acled.4$coefficients[control]/
                               lm.acled.4$se[control], dof = dof)
    }
    if(model[controlvar.index] == "lm.ucdp.4"){
      dof <- degrees_freedom(lm.ucdp.4, type = "resid")
      r2yxj.dx <- partial_r2(t_statistic = lm.ucdp.4$coefficients[control]/
                               lm.ucdp.4$se[control], dof = dof)
    }
    
    
# Use the t-value of control variable in the *treatment* regression
# to compute the partial R2 of control variable with the treatment
    if(model[controlvar.index] == "lm.acled.4"){
      if(control == "sqkm_ln"){
        formula_aux_reg <- paste(treatvar[controlvar.index], " ~ ", control, " + ",
                                 "policing_subgroups + mean_neigh_acled +",
                                 "status_egip + capdist_wmean_ex_ln + apop_ln +",
                                 "nightlight_corr_ex_ln + ",
                                 "vdem_egaldem + incidence_flag | year + cow")
      }
      if(control == "apop_ln"){
        formula_aux_reg <- paste(treatvar[controlvar.index], " ~ ", controlvar, " + ",
                                 "policing_subgroups + mean_neigh_acled +",
                                 "status_egip + capdist_wmean_ex_ln +",
                                 "sqkm_ln + nightlight_corr_ex_ln + ",
                                 "vdem_egaldem + incidence_flag | year + cow")
      }
    }
    
    if(model[controlvar.index] == "lm.ucdp.4"){
      if(control == "sqkm_ln"){
        formula_aux_reg <- paste(treatvar[controlvar.index], " ~ ", control, " + ",
                                 "policing_subgroups + mean_neigh_ucdp +",
                                 "status_egip + capdist_wmean_ex_ln + apop_ln +",
                                 "nightlight_corr_ex_ln + ",
                                 "vdem_egaldem + incidence_flag | year + cow")
      }
      if(control == "apop_ln"){
        formula_aux_reg <- paste(treatvar[controlvar.index], " ~ ", controlvar, " + ",
                                 "policing_subgroups + mean_neigh_ucdp +",
                                 "status_egip + capdist_wmean_ex_ln +",
                                 "sqkm_ln + nightlight_corr_ex_ln + ",
                                 "vdem_egaldem + incidence_flag | year + cow")
      }
    }
    
    
    aux_reg <- feols(as.formula(formula_aux_reg),
                     data = nc_tpi, vcov = "hetero")
    dof_treatment <- degrees_freedom(aux_reg, type = "resid")
    
# R2 of control variable
    r2dxj.x <- partial_r2(t_statistic = aux_reg$coefficients[control]/
                            aux_reg$se[control], 
                          dof = dof_treatment)
    

# Compute bounds on the strength of confounders manually
    
    bounds <- ovb_partial_r2_bound(r2dxj.x = r2dxj.x,
                                   r2yxj.dx = r2yxj.dx,
                                   kd = k,
                                   ky = k,
                                   bound_label = 
                                     paste(k, "x", paste(labels[controlvar.index])))
    
# Compute adjusted estimates manually
    if(model[controlvar.index] == "lm.acled.4"){
      if(type[controlvar.index] == "point"){
        bound.values <- adjusted_estimate(estimate = 
                                            lm.acled.4$coefficients[treatvar[controlvar.index]],
                                          se = lm.acled.4$se[treatvar[controlvar.index]],
                                          dof = dof,
                                          r2dz.x = bounds$r2dz.x,
                                          r2yz.dx = bounds$r2yz.dx)
      }
      if(type[controlvar.index] == "t"){
        bound.values <- adjusted_t(estimate = 
                                     lm.acled.4$coefficients[treatvar[controlvar.index]],
                                   se = lm.acled.4$se[treatvar[controlvar.index]],
                                   dof = dof,
                                   r2dz.x = bounds$r2dz.x,
                                   r2yz.dx = bounds$r2yz.dx)
      }
    }
    
    if(model[controlvar.index] == "lm.ucdp.4"){
      if(type[controlvar.index] == "point"){
        bound.values <- adjusted_estimate(estimate = 
                                            lm.ucdp.4$coefficients[treatvar[controlvar.index]],
                                          se = lm.ucdp.4$se[treatvar[controlvar.index]],
                                          dof = dof,
                                          r2dz.x = bounds$r2dz.x,
                                          r2yz.dx = bounds$r2yz.dx)
      }
      if(type[controlvar.index] == "t"){
        bound.values <- adjusted_t(estimate = 
                                     lm.ucdp.4$coefficients[treatvar[controlvar.index]],
                                   se = lm.ucdp.4$se[treatvar[controlvar.index]],
                                   dof = dof,
                                   r2dz.x = bounds$r2dz.x,
                                   r2yz.dx = bounds$r2yz.dx)
      }
    }
  }
  return(list(bounds, bound.values, dof))
}  


## ACLED point estimate:

out_settlement <- plot_sensitivity(labels = c("Settlement Area (log)"),
                                   controlvar = c("sqkm_ln"),
                                   treatvar = c("policing_capacity"),
                                   model = c("lm.acled.4"),
                                   type = c("point"))

out_popsize <- plot_sensitivity(labels = c("Population Size (log)"),
                                controlvar = c("apop_ln"),
                                treatvar = c("policing_capacity"),
                                model = c("lm.acled.4"),
                                type = c("point"))

# Plot contours and bounds
pdf(file= "figures/sense_acled_cap.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = lm.acled.4$coefficients["policing_capacity"],
                 se = lm.acled.4$se["policing_capacity"],
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)
add_bound_to_contour(out_settlement[[1]], bound_value = out_settlement[[2]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], bound_value = out_popsize[[2]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, "(-0.013)", cex=1.5)
text(0.07, 0.044, "3 x Population Size (log)", cex=1.5)
text(0.07, 0.034, "(-0.009)", cex=1.5)
text(0.12, 0.08, "3 x Settlement Area (log)", cex=1.5)
text(0.12, 0.07, "(0.005)", cex=1.5)

dev.off()


## ACLED t-value:
out_settlement <- plot_sensitivity(labels = c("Settlement Area (log)"),
                                   controlvar = c("sqkm_ln"),
                                   treatvar = c("policing_capacity"),
                                   model = c("lm.acled.4"),
                                   type = c("t"))

out_popsize <- plot_sensitivity(labels = c("Population Size (log)"),
                                controlvar = c("apop_ln"),
                                treatvar = c("policing_capacity"),
                                model = c("lm.acled.4"),
                                type = c("t"))

pdf(file= "figures/sense_acled_cap_t.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = lm.acled.4$coefficients["policing_capacity"],
                 se = lm.acled.4$se["policing_capacity"],
                 sensitivity.of = "t-value",
                 t.threshold = -1.96,
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)
add_bound_to_contour(out_settlement[[1]], bound_value = out_settlement[[2]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], bound_value = out_popsize[[2]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

t.unadj <- lm.acled.4$coefficients["policing_capacity"]/lm.acled.4$se["policing_capacity"]
t.unadj <- signif(t.unadj, 3)
t.unadj <- paste("(", t.unadj, ")", sep = "")

t.settlement <- signif(out_settlement[[2]], 1)
t.settlement <- paste("(", t.settlement, ")", sep = "")

t.pop <- signif(out_popsize[[2]], 2)
t.pop <- paste("(", t.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, t.unadj, cex=1.5)
text(0.07, 0.042, "3 x Population Size (log)", cex=1.5)
text(0.07, 0.032, t.pop, cex=1.5)
text(0.11, 0.07, "3 x Settlement Area (log)", cex=1.5)
text(0.11, 0.06, t.settlement, cex=1.5)
dev.off() 


### UCDP Point Estimate

out_settlement <- plot_sensitivity(labels = c("Settlement Area (log)"),
                                   controlvar = c("sqkm_ln"),
                                   treatvar = c("policing_capacity"),
                                   model = c("lm.ucdp.4"),
                                   type = c("point"),
                                   k = 6)

out_popsize <- plot_sensitivity(labels = c("Population Size (log)"),
                                controlvar = c("apop_ln"),
                                treatvar = c("policing_capacity"),
                                model = c("lm.ucdp.4"),
                                type = c("point"),
                                k = 6)

pdf(file= "figures/sense_ucdp_cap.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = lm.ucdp.4$coefficients["policing_capacity"],
                 se = lm.ucdp.4$se["policing_capacity"],
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)
add_bound_to_contour(out_settlement[[1]], bound_value = signif(out_settlement[[2]], digits = 1),
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], bound_value = signif(out_popsize[[2]], digits = 1),
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- lm.ucdp.4$coefficients["policing_capacity"]
est.unadj <- signif(est.unadj, 1)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 1)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 1)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, est.unadj , cex=1.5)
text(0.08, 0.03, "6 x Population Size (log)", cex=1.5)
text(0.08, 0.02, est.pop, cex=1.5)
text(0.20, 0.015, "6 x Settlement Area (log) ", cex=1.5)
text(0.20, 0.005, est.settlement , cex=1.5)
dev.off() 


## UCDP t-value

out_settlement <- plot_sensitivity(labels = c("Settlement Area (log)"),
                                   controlvar = c("sqkm_ln"),
                                   treatvar = c("policing_capacity"),
                                   model = c("lm.ucdp.4"),
                                   type = c("t"),
                                   k = 6)

out_popsize <- plot_sensitivity(labels = c("Population Size (log)"),
                                controlvar = c("apop_ln"),
                                treatvar = c("policing_capacity"),
                                model = c("lm.ucdp.4"),
                                type = c("t"),
                                k = 6)

pdf(file= "figures/sense_ucdp_cap_t.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = lm.ucdp.4$coefficients["policing_capacity"],
                 se = lm.ucdp.4$se["policing_capacity"],
                 sensitivity.of = "t-value",
                 t.threshold = -1.96,
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)
add_bound_to_contour(out_settlement[[1]], bound_value = signif(out_settlement[[2]], digits = 1),
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], bound_value = signif(out_popsize[[2]], digits = 1),
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- lm.ucdp.4$coefficients["policing_capacity"]/lm.ucdp.4$se["policing_capacity"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, est.unadj, cex=1.5)
text(0.08, 0.03, "6 x Population Size (log)", cex=1.5)
text(0.08, 0.02, est.pop, cex=1.5)
text(0.22, 0.05, "6 x Settlement Area (log) ", cex=1.5)
text(0.22, 0.04, est.settlement, cex=1.5)
dev.off() 


## > Figure A7: Sensitivity Direct Effect of Policing Capacity ----

# Relevant models are specified above:
  # med.lm.acled.4
  # med.lm.ucdp.4


# preparatory function for sensemakr with direct effect models:
plot_sensitivity_direct <- function(controlvar, treatvar,
                                    model, labels, type, k = 3){
  
  for(control in controlvar){
    
    controlvar.index <- which(controlvar == control)
    
# Use the t statistic of control variable in the outcome regression
# to compute the partial R2 of control variable with the outcome.
    if(model[controlvar.index] == "med.lm.acled.4"){
      dof <- degrees_freedom(med.lm.acled.4, type = "resid")
      r2yxj.dx <- partial_r2(t_statistic = med.lm.acled.4$coefficients[control]/
                               med.lm.acled.4$se[control], dof = dof)
    }
    if(model[controlvar.index] == "med.lm.ucdp.4"){
      dof <- degrees_freedom(med.lm.ucdp.4, type = "resid")
      r2yxj.dx <- partial_r2(t_statistic = med.lm.ucdp.4$coefficients[control]/
                               med.lm.ucdp.4$se[control], dof = dof)
    }
    
    
# Use the t-value of control variable in the *treatment* regression
# to compute the partial R2 of control variable with the treatment
    if(model[controlvar.index] == "med.lm.acled.4"){
      if(control == "sqkm_ln"){
        formula_aux_reg <- paste(treatvar[controlvar.index], " ~ ", control, " + ",
                                 "policing_subgroups + mean_neigh_acled +",
                                 "policing_conflict + ",
                                 "status_egip + capdist_wmean_ex_ln + apop_ln +",
                                 "nightlight_corr_ex_ln + ",
                                 "vdem_egaldem + incidence_flag | year + cow")
      }
      if(control == "apop_ln"){
        formula_aux_reg <- paste(treatvar[controlvar.index], " ~ ", controlvar, " + ",
                                 "policing_subgroups + mean_neigh_acled +",
                                 "status_egip + capdist_wmean_ex_ln +",
                                 "policing_conflict + ",
                                 "sqkm_ln + nightlight_corr_ex_ln + ",
                                 "vdem_egaldem + incidence_flag | year + cow")
      }
    }
    
    if(model[controlvar.index] == "med.lm.ucdp.4"){
      if(control == "sqkm_ln"){
        formula_aux_reg <- paste(treatvar[controlvar.index], " ~ ", control, " + ",
                                 "policing_subgroups + mean_neigh_ucdp +",
                                 "policing_conflict + ",
                                 "status_egip + capdist_wmean_ex_ln + apop_ln +",
                                 "nightlight_corr_ex_ln + ",
                                 "vdem_egaldem + incidence_flag | year + cow")
      }
      if(control == "apop_ln"){
        formula_aux_reg <- paste(treatvar[controlvar.index], " ~ ", controlvar, " + ",
                                 "policing_subgroups + mean_neigh_ucdp +",
                                 "policing_conflict + ",
                                 "status_egip + capdist_wmean_ex_ln +",
                                 "sqkm_ln + nightlight_corr_ex_ln + ",
                                 "vdem_egaldem + incidence_flag | year + cow")
      }
    }
    
    
    aux_reg <- feols(as.formula(formula_aux_reg),
                     data = nc_tpi, vcov = "hetero")
    dof_treatment <- degrees_freedom(aux_reg, type = "resid")
    
# R2 of control variable
    r2dxj.x <- partial_r2(t_statistic = aux_reg$coefficients[control]/
                            aux_reg$se[control], 
                          dof = dof_treatment)
    
# Compute bounds on the strength of confounders manually
    
    bounds <- ovb_partial_r2_bound(r2dxj.x = r2dxj.x,
                                   r2yxj.dx = r2yxj.dx,
                                   kd = k,
                                   ky = k,
                                   bound_label = 
                                     paste(k, "x", paste(labels[controlvar.index])))
    
# Compute adjusted estimates manually
    if(model[controlvar.index] == "med.lm.acled.4"){
      if(type[controlvar.index] == "point"){
        bound.values <- adjusted_estimate(estimate = 
                                            med.lm.acled.4$coefficients[treatvar[controlvar.index]],
                                          se = med.lm.acled.4$se[treatvar[controlvar.index]],
                                          dof = dof,
                                          r2dz.x = bounds$r2dz.x,
                                          r2yz.dx = bounds$r2yz.dx)
      }
      if(type[controlvar.index] == "t"){
        bound.values <- adjusted_t(estimate = 
                                     med.lm.acled.4$coefficients[treatvar[controlvar.index]],
                                   se = med.lm.acled.4$se[treatvar[controlvar.index]],
                                   dof = dof,
                                   r2dz.x = bounds$r2dz.x,
                                   r2yz.dx = bounds$r2yz.dx)
      }
    }
    
    if(model[controlvar.index] == "med.lm.ucdp.4"){
      if(type[controlvar.index] == "point"){
        bound.values <- adjusted_estimate(estimate = 
                                            med.lm.ucdp.4$coefficients[treatvar[controlvar.index]],
                                          se = med.lm.ucdp.4$se[treatvar[controlvar.index]],
                                          dof = dof,
                                          r2dz.x = bounds$r2dz.x,
                                          r2yz.dx = bounds$r2yz.dx)
      }
      if(type[controlvar.index] == "t"){
        bound.values <- adjusted_t(estimate = 
                                     med.lm.ucdp.4$coefficients[treatvar[controlvar.index]],
                                   se = med.lm.ucdp.4$se[treatvar[controlvar.index]],
                                   dof = dof,
                                   r2dz.x = bounds$r2dz.x,
                                   r2yz.dx = bounds$r2yz.dx)
      }
    }
  }
  return(list(bounds, bound.values, dof))
}



## ACLED point estimate

out_settlement <- plot_sensitivity_direct(labels = "Settlement Area (log)",
                                          controlvar = "sqkm_ln",
                                          treatvar = "policing_capacity",
                                          model = "med.lm.acled.4",
                                          type = "point",
                                          k = 3)

out_popsize <- plot_sensitivity_direct(labels = c("Population Size (log)"),
                                       controlvar = c("apop_ln"),
                                       treatvar = c("policing_capacity"),
                                       model = c("med.lm.acled.4"),
                                       type = c("point"),
                                       k = 3)

pdf(file= "figures/sense_med_acled_cap.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = med.lm.acled.4$coefficients["policing_capacity"],
                 se = med.lm.acled.4$se["policing_capacity"],
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)
add_bound_to_contour(out_settlement[[1]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], 
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- med.lm.acled.4$coefficients["policing_capacity"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, est.unadj, cex=1.5)
text(0.06, 0.06, "3 x Settlement Area (log)", cex=1.5)
text(0.06, 0.05, est.settlement , cex=1.5)
text(0.10, 0.0275, "3 x Population Size (log)", cex=1.5)
text(0.10, 0.0175, est.pop, cex=1.5)

dev.off() 


## ACLED t-value

out_settlement <- plot_sensitivity_direct(labels = "Settlement Area (log)",
                                          controlvar = "sqkm_ln",
                                          treatvar = "policing_capacity",
                                          model = "med.lm.acled.4",
                                          type = "t",
                                          k = 3)

out_popsize <- plot_sensitivity_direct(labels = c("Population Size (log)"),
                                       controlvar = c("apop_ln"),
                                       treatvar = c("policing_capacity"),
                                       model = c("med.lm.acled.4"),
                                       type = c("t"),
                                       k = 3)

pdf(file= "figures/sense_med_acled_cap_t.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = med.lm.acled.4$coefficients["policing_capacity"],
                 se = med.lm.acled.4$se["policing_capacity"],
                 sensitivity.of = "t-value",
                 t.threshold = -1.96,
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)
add_bound_to_contour(out_settlement[[1]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], 
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- med.lm.acled.4$coefficients["policing_capacity"]/med.lm.acled.4$se["policing_capacity"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, est.unadj, cex=1.5)
text(0.07, 0.05, "3 x Settlement Area (log)", cex=1.5)
text(0.07, 0.04, est.settlement , cex=1.5)
text(0.09, 0.03, "3 x Population Size (log)", cex=1.5)
text(0.09, 0.02, est.pop, cex=1.5)
dev.off() 


## UCDP point estimate

out_settlement <- plot_sensitivity_direct(labels = "Settlement Area (log)",
                                          controlvar = "sqkm_ln",
                                          treatvar = "policing_capacity",
                                          model = "med.lm.ucdp.4",
                                          type = "point",
                                          k = 3)

out_popsize <- plot_sensitivity_direct(labels = c("Population Size (log)"),
                                       controlvar = c("apop_ln"),
                                       treatvar = c("policing_capacity"),
                                       model = c("med.lm.ucdp.4"),
                                       type = c("point"),
                                       k = 3)

pdf(file= "figures/sense_med_ucdp_cap.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = med.lm.ucdp.4$coefficients["policing_capacity"],
                 se = med.lm.ucdp.4$se["policing_capacity"],
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)
add_bound_to_contour(out_settlement[[1]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], 
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- med.lm.ucdp.4$coefficients["policing_capacity"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, est.unadj , cex=1.5)
text(0.15, 0.025, "3 x Settlement Area (log) / 3 x Population Size (log)", cex=1.5)
text(0.15, 0.015, "(-0.008) / (-0.007)", cex=1.5)
dev.off() 


## UCDP t-value
out_settlement <- plot_sensitivity_direct(labels = "Settlement Area (log)",
                                          controlvar = "sqkm_ln",
                                          treatvar = "policing_capacity",
                                          model = "med.lm.ucdp.4",
                                          type = "t",
                                          k = 3)

out_popsize <- plot_sensitivity_direct(labels = c("Population Size (log)"),
                                       controlvar = c("apop_ln"),
                                       treatvar = c("policing_capacity"),
                                       model = c("med.lm.ucdp.4"),
                                       type = c("t"),
                                       k = 3)

pdf(file= "figures/sense_med_ucdp_cap_t.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = med.lm.ucdp.4$coefficients["policing_capacity"],
                 se = med.lm.ucdp.4$se["policing_capacity"],
                 sensitivity.of = "t-value",
                 t.threshold = -1.96,
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)
add_bound_to_contour(out_settlement[[1]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], 
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- med.lm.ucdp.4$coefficients["policing_capacity"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, "(-4.87)", cex=1.5)
text(0.17, 0.025, "3 x Settlement Area (log) / 3 x Population Size (log)", cex=1.5)
text(0.17, 0.015, "(-2.361) / (-1.996)", cex=1.5)
dev.off() 




## > Figure A8: Sensitivity of Subgroup Count with Policing ----

#Notes: 
  # Relevant models are specified above:
    # lm.acled.4
    # lm.ucdp.4
  #Run plot_sensitivity function above (Figure A6) for the code below


# ACLED point estimate

out_settlement <- plot_sensitivity(labels = "Settlement Area (log)",
                                   controlvar = "sqkm_ln",
                                   treatvar = "policing_subgroups",
                                   model = "lm.acled.4",
                                   type = "point",
                                   k = 3)

out_popsize <- plot_sensitivity(labels = c("Population Size (log)"),
                                controlvar = c("apop_ln"),
                                treatvar = c("policing_subgroups"),
                                model = c("lm.acled.4"),
                                type = c("point"),
                                k = 3)

pdf(file= "figures/sense_acled_horiz.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = lm.acled.4$coefficients["policing_subgroups"],
                 se = lm.acled.4$se["policing_subgroups"],
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)

add_bound_to_contour(out_settlement[[1]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], 
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- lm.acled.4$coefficients["policing_subgroups"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, est.unadj, cex=1.5)
text(0.135, 0.03, "3 x Population Size (log)", cex=1.5)
text(0.135, 0.02, est.pop, cex=1.5)
text(0.11, 0.065, "3 x Settlement Area (log)", cex=1.5)
text(0.11, 0.055, est.settlement, cex=1.5)
dev.off() 


# ACLED t-value

out_settlement <- plot_sensitivity(labels = "Settlement Area (log)",
                                   controlvar = "sqkm_ln",
                                   treatvar = "policing_subgroups",
                                   model = "lm.acled.4",
                                   type = "t",
                                   k = 3)

out_popsize <- plot_sensitivity(labels = c("Population Size (log)"),
                                controlvar = c("apop_ln"),
                                treatvar = c("policing_subgroups"),
                                model = c("lm.acled.4"),
                                type = c("t"),
                                k = 3)

pdf(file= "figures/sense_acled_horiz_t.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = lm.acled.4$coefficients["policing_subgroups"],
                 se = lm.acled.4$se["policing_subgroups"],
                 sensitivity.of = "t-value",
                 t.threshold = 1.96,
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)

add_bound_to_contour(out_settlement[[1]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], 
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- lm.acled.4$coefficients["policing_subgroups"]/lm.acled.4$se["policing_subgroups"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0.005, "Unadjusted", cex=1.5)
text(0.032, -0.005, est.unadj, cex=1.5)
text(0.135, 0.03, "3 x Population Size (log)", cex=1.5)
text(0.135, 0.02, est.pop, cex=1.5)
text(0.12, 0.065, "3 x Settlement Area (log)", cex=1.5)
text(0.12, 0.055, est.settlement, cex=1.5)
dev.off() 



# UCDP point estimate

out_settlement <- plot_sensitivity(labels = "Settlement Area (log)",
                                   controlvar = "sqkm_ln",
                                   treatvar = "policing_subgroups",
                                   model = "lm.ucdp.4",
                                   type = "point",
                                   k = 3)

out_popsize <- plot_sensitivity(labels = c("Population Size (log)"),
                                controlvar = c("apop_ln"),
                                treatvar = c("policing_subgroups"),
                                model = c("lm.ucdp.4"),
                                type = c("point"),
                                k = 3)

pdf(file= "figures/sense_ucdp_horiz.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = lm.ucdp.4$coefficients["policing_subgroups"],
                 se = lm.ucdp.4$se["policing_subgroups"],
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)

add_bound_to_contour(out_settlement[[1]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], 
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- lm.ucdp.4$coefficients["policing_subgroups"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0, "Unadjusted", cex=1.5)
text(0.032, -0.01, est.unadj, cex=1.5)
text(0.175, 0.0225, "6 x Population Size (log)", cex=1.5)
text(0.175, 0.0125, est.pop, cex=1.5)
text(0.045, 0.025, "6 x Settlement Area (log)", cex=1.5)
text(0.045, 0.015, est.settlement, cex=1.5)
dev.off() 


# UCDP t-value

out_settlement <- plot_sensitivity(labels = "Settlement Area (log)",
                                   controlvar = "sqkm_ln",
                                   treatvar = "policing_subgroups",
                                   model = "lm.ucdp.4",
                                   type = "t",
                                   k = 3)

out_popsize <- plot_sensitivity(labels = c("Population Size (log)"),
                                controlvar = c("apop_ln"),
                                treatvar = c("policing_subgroups"),
                                model = c("lm.ucdp.4"),
                                type = c("t"),
                                k = 3)

pdf(file= "figures/sense_ucdp_horiz_t.pdf",  
    height=11, width=12)

ovb_contour_plot(estimate = lm.ucdp.4$coefficients["policing_subgroups"],
                 se = lm.ucdp.4$se["policing_subgroups"],
                 sensitivity.of = "t-value",
                 t.threshold = 1.96,
                 dof = out_settlement[[3]],
                 label.text	= F,
                 cex.lab=1.5,
                 labcex = 1.5,
                 cex.axis = 1.5,
                 cex.label.text = 0.001)

add_bound_to_contour(out_settlement[[1]],
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

add_bound_to_contour(out_popsize[[1]], 
                     label.text	= F,
                     cex.lab=1.5,
                     labcex = 1.5,
                     cex.axis = 1.5,
                     cex.label.text = 0.001)

est.unadj <- lm.ucdp.4$coefficients["policing_subgroups"]/lm.ucdp.4$se["policing_subgroups"]
est.unadj <- signif(est.unadj, 3)
est.unadj <- paste("(", est.unadj, ")", sep = "")

est.settlement <- signif(out_settlement[[2]], 2)
est.settlement <- paste("(", est.settlement, ")", sep = "")

est.pop <- signif(out_popsize[[2]], 2)
est.pop <- paste("(", est.pop, ")", sep = "")

text(0.032, 0, "Unadjusted", cex=1.5)
text(0.032, -0.01, est.unadj, cex=1.5)
text(0.18, 0.02, "3 x Population Size (log)", cex=1.5)
text(0.18, 0.01, est.pop, cex=1.5)
text(0.045, 0.025, "3 x Settlement Area (log)", cex=1.5)
text(0.045, 0.015, est.settlement, cex=1.5)
dev.off() 




## > Tables A6 & A7 ----

# ACLED Supergroup Policing
fileConn<-file(file.path("tables", paste0("sense_acled_cap", ".tex")))
writeLines(ovb_minimal_reporting(
  sensemakr:::sensemakr.numeric(
    estimate = lm.acled.4$coefficients["policing_capacity"],
    se = lm.acled.4$se["policing_capacity"],
    dof = degrees_freedom(lm.acled.4, type = "resid")),
  format = "latex"), 
  fileConn)
close(fileConn)

# UCDP Supergroup Policing
fileConn<-file(file.path("tables", paste0("sense_ucdp_cap", ".tex")))
writeLines(ovb_minimal_reporting(
  sensemakr:::sensemakr.numeric(
    estimate = lm.ucdp.4$coefficients["policing_capacity"],
    se = lm.ucdp.4$se["policing_capacity"],
    dof = degrees_freedom(lm.ucdp.4, type = "resid")),
  format = "latex"), 
  fileConn)
close(fileConn)


# ACLED Supergroup Policing Direct Effect
fileConn<-file(file.path("tables", paste0("sense_med_acled_cap", ".tex")))
writeLines(ovb_minimal_reporting(
  sensemakr:::sensemakr.numeric(
    estimate = med.lm.acled.4$coefficients["policing_capacity"],
    se = med.lm.acled.4$se["policing_capacity"],
    dof = degrees_freedom(med.lm.acled.4, type = "resid")),
  format = "latex"), 
  fileConn)
close(fileConn)

# UCDP Supergroup Policing Direct Effect
fileConn<-file(file.path("tables", paste0("sense_med_ucdp_cap", ".tex")))
writeLines(ovb_minimal_reporting(
  sensemakr:::sensemakr.numeric(
    estimate = med.lm.ucdp.4$coefficients["policing_capacity"],
    se = med.lm.ucdp.4$se["policing_capacity"],
    dof = degrees_freedom(med.lm.ucdp.4, type = "resid")),
  format = "latex"), 
  fileConn)
close(fileConn)


# ACLED No. of Subgroups w/ Policing
fileConn<-file(file.path("tables", paste0("sense_acled_horiz", ".tex")))
writeLines(ovb_minimal_reporting(
  sensemakr:::sensemakr.numeric(
    estimate = lm.acled.4$coefficients["policing_subgroups"],
    se = lm.acled.4$se["policing_subgroups"],
    dof = degrees_freedom(lm.acled.4, type = "resid")),
  format = "latex"), 
  fileConn)
close(fileConn)

# UCDP No. of Subgroups w/ Policing
fileConn<-file(file.path("tables", paste0("sense_ucdp_horiz", ".tex")))
writeLines(ovb_minimal_reporting(
  sensemakr:::sensemakr.numeric(
    estimate = lm.ucdp.4$coefficients["policing_subgroups"],
    se = lm.ucdp.4$se["policing_subgroups"],
    dof = degrees_freedom(lm.ucdp.4, type = "resid")),
  format = "latex"), 
  fileConn)
close(fileConn)






####### Figure A3: Robustness ----

### > Mixed Effect ----

# lmer functions for main effect
formulae.me.1 <- sapply(outcomes_lags,
                        function(x) reformulate(c(treat.1, x[2], baseline.me), response = x[1]))
formulae.me.2 <- sapply(outcomes_lags, 
                        function(x) reformulate(c(treat.2, x[2], baseline.me), response = x[1]))
formulae.me.3 <- sapply(outcomes_lags, 
                        function(x) reformulate(c(treat, x[2], baseline.me), response = x[1]))
formulae.me.4 <- sapply(outcomes_lags, 
                        function(x) reformulate(c(treat, x[2], baseline.me, controls.me), response = x[1]))

# UCDP
model.ls <- list(
  lmer.ucdp.1 <- lmer(formulae.me.1[[2]], 
                      data = nc_tpi),
  lmer.ucdp.2 <- lmer(formulae.me.2[[2]], 
                      data = nc_tpi),
  lmer.ucdp.3 <- lmer(formulae.me.3[[2]], 
                      data = nc_tpi),
  lmer.ucdp.4 <- lmer(formulae.me.4[[2]], 
                      data = nc_tpi)
)

# ACLED
model.ls <- list(
  lmer.acled.1 <- lmer(formulae.me.1[[1]], 
                       data = nc_tpi),
  lmer.acled.2 <- lmer(formulae.me.2[[1]], 
                       data = nc_tpi),
  lmer.acled.3 <- lmer(formulae.me.3[[1]], 
                       data = nc_tpi),
  lmer.acled.4 <- lmer(formulae.me.4[[1]], 
                       data = nc_tpi)
)


#Extract estimates
estimates.lmer_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Mixed Effect","Mixed Effect","Mixed Effect","Mixed Effect","Mixed Effect","Mixed Effect"), 
  estimate = c(coef(lmer.ucdp.2)$cow[1,"policing_subgroups"],
               coef(lmer.ucdp.3)$cow[1,"policing_subgroups"],
               coef(lmer.ucdp.4)$cow[1,"policing_subgroups"],
               coef(lmer.ucdp.1)$cow[1,"policing_capacity"],
               coef(lmer.ucdp.3)$cow[1,"policing_capacity"],
               coef(lmer.ucdp.4)$cow[1,"policing_capacity"]),
  ci.lower = c(confint(lmer.ucdp.2)["policing_subgroups", 1],
               confint(lmer.ucdp.3)["policing_subgroups", 1],
               confint(lmer.ucdp.4)["policing_subgroups", 1],
               confint(lmer.ucdp.1)["policing_capacity", 1],
               confint(lmer.ucdp.3)["policing_capacity", 1],
               confint(lmer.ucdp.4)["policing_capacity", 1]),
  ci.upper = c(confint(lmer.ucdp.2)["policing_subgroups", 2],
               confint(lmer.ucdp.3)["policing_subgroups", 2],
               confint(lmer.ucdp.4)["policing_subgroups", 2],
               confint(lmer.ucdp.1)["policing_capacity", 2],
               confint(lmer.ucdp.3)["policing_capacity", 2],
               confint(lmer.ucdp.4)["policing_capacity", 2])   )  


estimates.lmer_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Mixed Effect","Mixed Effect","Mixed Effect","Mixed Effect","Mixed Effect","Mixed Effect"), 
  estimate = c(coef(lmer.acled.2)$cow[1,"policing_subgroups"],
               coef(lmer.acled.3)$cow[1,"policing_subgroups"],
               coef(lmer.acled.4)$cow[1,"policing_subgroups"],
               coef(lmer.acled.1)$cow[1,"policing_capacity"],
               coef(lmer.acled.3)$cow[1,"policing_capacity"],
               coef(lmer.acled.4)$cow[1,"policing_capacity"]),
  ci.lower = c(confint(lmer.acled.2)["policing_subgroups", 1],
               confint(lmer.acled.3)["policing_subgroups", 1],
               confint(lmer.acled.4)["policing_subgroups", 1],
               confint(lmer.acled.1)["policing_capacity", 1],
               confint(lmer.acled.3)["policing_capacity", 1],
               confint(lmer.acled.4)["policing_capacity", 1]),
  ci.upper = c(confint(lmer.acled.2)["policing_subgroups", 2],
               confint(lmer.acled.3)["policing_subgroups", 2],
               confint(lmer.acled.4)["policing_subgroups", 2],
               confint(lmer.acled.1)["policing_capacity", 2],
               confint(lmer.acled.3)["policing_capacity", 2],
               confint(lmer.acled.4)["policing_capacity", 2])   )  


estimates.lmer<- bind_rows(
  list(estimates.lmer_ucdp,
       estimates.lmer_acled)
)

estimates.lmer$model_factor <- factor(estimates.lmer$model, 
                                      levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.lmer, file="estimates/estimates.lmer.RData" )

### > All Groups ----


# ACLED
model.ls <- list(
  all_grps.lm.acled.1 <- feols(formulae$twfe$twfe.1[[10]],
                               data = nc_rep, vcov = "hetero"),
  all_grps.lm.acled.2 <- feols(formulae$twfe$twfe.2[[10]],
                               data = nc_rep, vcov = "hetero"),
  all_grps.lm.acled.3 <- feols(formulae$twfe$twfe.3[[10]],
                               data = nc_rep, vcov = "hetero"),
  all_grps.lm.acled.4 <- feols(formulae$twfe$twfe.4[[10]],
                               data = nc_rep, vcov = "hetero")
)


# UCDP
model.ls <- list(
  all_grps.lm.ucdp.1 <- feols(formulae$twfe$twfe.1[[11]],
                              data = nc_rep, vcov = "hetero"),
  all_grps.lm.ucdp.2 <- feols(formulae$twfe$twfe.2[[11]],
                              data = nc_rep, vcov = "hetero"),
  all_grps.lm.ucdp.3 <- feols(formulae$twfe$twfe.3[[11]],
                              data = nc_rep, vcov = "hetero"),
  all_grps.lm.ucdp.4 <- feols(formulae$twfe$twfe.4[[11]],
                              data = nc_rep, vcov = "hetero")
)


#Extract estimates

estimates.all_grps_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("All Groups", "All Groups","All Groups","All Groups","All Groups","All Groups"),
  estimate = c(coef(all_grps.lm.acled.2)["policing_subgroups"],
               coef(all_grps.lm.acled.3)["policing_subgroups"],
               coef(all_grps.lm.acled.4)["policing_subgroups"],
               coef(all_grps.lm.acled.1)["policing_capacity"],
               coef(all_grps.lm.acled.3)["policing_capacity"],
               coef(all_grps.lm.acled.4)["policing_capacity"]),
  ci.lower = c(confint(all_grps.lm.acled.2)["policing_subgroups", 1],
               confint(all_grps.lm.acled.3)["policing_subgroups", 1],
               confint(all_grps.lm.acled.4)["policing_subgroups", 1],
               confint(all_grps.lm.acled.1)["policing_capacity", 1],
               confint(all_grps.lm.acled.3)["policing_capacity", 1],
               confint(all_grps.lm.acled.4)["policing_capacity", 1]),
  ci.upper = c(confint(all_grps.lm.acled.2)["policing_subgroups", 2],
               confint(all_grps.lm.acled.3)["policing_subgroups", 2],
               confint(all_grps.lm.acled.4)["policing_subgroups", 2],
               confint(all_grps.lm.acled.1)["policing_capacity", 2],
               confint(all_grps.lm.acled.3)["policing_capacity", 2],
               confint(all_grps.lm.acled.4)["policing_capacity", 2])   )   


estimates.all_grps_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("All Groups", "All Groups","All Groups","All Groups","All Groups","All Groups"),
  estimate = c(coef(all_grps.lm.ucdp.2)["policing_subgroups"],
               coef(all_grps.lm.ucdp.3)["policing_subgroups"],
               coef(all_grps.lm.ucdp.4)["policing_subgroups"],
               coef(all_grps.lm.ucdp.1)["policing_capacity"],
               coef(all_grps.lm.ucdp.3)["policing_capacity"],
               coef(all_grps.lm.ucdp.4)["policing_capacity"]),
  ci.lower = c(confint(all_grps.lm.ucdp.2)["policing_subgroups", 1],
               confint(all_grps.lm.ucdp.3)["policing_subgroups", 1],
               confint(all_grps.lm.ucdp.4)["policing_subgroups", 1],
               confint(all_grps.lm.ucdp.1)["policing_capacity", 1],
               confint(all_grps.lm.ucdp.3)["policing_capacity", 1],
               confint(all_grps.lm.ucdp.4)["policing_capacity", 1]),
  ci.upper = c(confint(all_grps.lm.ucdp.2)["policing_subgroups", 2],
               confint(all_grps.lm.ucdp.3)["policing_subgroups", 2],
               confint(all_grps.lm.ucdp.4)["policing_subgroups", 2],
               confint(all_grps.lm.ucdp.1)["policing_capacity", 2],
               confint(all_grps.lm.ucdp.3)["policing_capacity", 2],
               confint(all_grps.lm.ucdp.4)["policing_capacity", 2])   )  



## Interaction

# ACLED
model.ls <- list(
  all_grps.inter.lmer.acled.1 <- lmer(formulae.me.inter.1[10, 1][[1]], 
                                      data = nc_rep),
  all_grps.inter.lmer.acled.2 <- lmer(formulae.me.inter.2[10, 1][[1]], 
                                      data = nc_rep),
  all_grps.inter.lmer.acled.3 <- lmer(formulae.me.inter.1[10, 2][[1]], 
                                      data = nc_rep),
  all_grps.inter.lmer.acled.4 <- lmer(formulae.me.inter.2[10, 2][[1]], 
                                      data = nc_rep)
)

# UCDP
model.ls <- list(
  all_grps.inter.lmer.ucdp.1 <- lmer(formulae.me.inter.1[11, 1][[1]], 
                                     data = nc_rep),
  all_grps.inter.lmer.ucdp.2 <- lmer(formulae.me.inter.2[11, 1][[1]], 
                                     data = nc_rep),
  all_grps.inter.lmer.ucdp.3 <- lmer(formulae.me.inter.1[11, 2][[1]], 
                                     data = nc_rep),
  all_grps.inter.lmer.ucdp.4 <- lmer(formulae.me.inter.2[11, 2][[1]], 
                                     data = nc_rep)
)



#Extract estimates
estimates.inter.all_grps <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("All Groups", "All Groups","All Groups","All Groups","All Groups","All Groups","All Groups","All Groups"),
  estimate = c(coef(all_grps.inter.lmer.acled.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(all_grps.inter.lmer.acled.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(all_grps.inter.lmer.acled.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(all_grps.inter.lmer.acled.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(all_grps.inter.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(all_grps.inter.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(all_grps.inter.lmer.ucdp.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(all_grps.inter.lmer.ucdp.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(all_grps.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(all_grps.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(all_grps.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(all_grps.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(all_grps.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(all_grps.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(all_grps.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(all_grps.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(all_grps.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(all_grps.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(all_grps.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(all_grps.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(all_grps.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(all_grps.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(all_grps.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(all_grps.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )

estimates.all_grps <- bind_rows(
  list(estimates.all_grps_ucdp,
       estimates.all_grps_acled,
       estimates.inter.all_grps)
)


estimates.all_grps$model_factor <- factor(estimates.all_grps$model, 
                                          levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.all_grps, file="estimates/estimates.all_grps.RData" )



### > Groups Outside Sub-Saharan Africa ----


# ACLED
model.ls <- list(
  nonAfr.lm.acled.1 <- feols(formulae$twfe$twfe.1[[1]],
                             data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.lm.acled.2 <- feols(formulae$twfe$twfe.2[[1]],
                             data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.lm.acled.3 <- feols(formulae$twfe$twfe.3[[1]],
                             data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.lm.acled.4 <- feols(formulae$twfe$twfe.4[[1]],
                             data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",])
)


# UCDP
model.ls <- list(
  nonAfr.lm.ucdp.1 <- feols(formulae$twfe$twfe.1[[2]],
                            data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.lm.ucdp.2 <- feols(formulae$twfe$twfe.2[[2]],
                            data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.lm.ucdp.3 <- feols(formulae$twfe$twfe.3[[2]],
                            data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.lm.ucdp.4 <- feols(formulae$twfe$twfe.4[[2]],
                            data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",])
)



#Extract estimates

estimates.nonAfr_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Groups outside \n Sub-Saharan Africa", "Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa"),
  estimate = c(coef(nonAfr.lm.acled.2)["policing_subgroups"],
               coef(nonAfr.lm.acled.3)["policing_subgroups"],
               coef(nonAfr.lm.acled.4)["policing_subgroups"],
               coef(nonAfr.lm.acled.1)["policing_capacity"],
               coef(nonAfr.lm.acled.3)["policing_capacity"],
               coef(nonAfr.lm.acled.4)["policing_capacity"]),
  ci.lower = c(confint(nonAfr.lm.acled.2)["policing_subgroups", 1],
               confint(nonAfr.lm.acled.3)["policing_subgroups", 1],
               confint(nonAfr.lm.acled.4)["policing_subgroups", 1],
               confint(nonAfr.lm.acled.1)["policing_capacity", 1],
               confint(nonAfr.lm.acled.3)["policing_capacity", 1],
               confint(nonAfr.lm.acled.4)["policing_capacity", 1]),
  ci.upper = c(confint(nonAfr.lm.acled.2)["policing_subgroups", 2],
               confint(nonAfr.lm.acled.3)["policing_subgroups", 2],
               confint(nonAfr.lm.acled.4)["policing_subgroups", 2],
               confint(nonAfr.lm.acled.1)["policing_capacity", 2],
               confint(nonAfr.lm.acled.3)["policing_capacity", 2],
               confint(nonAfr.lm.acled.4)["policing_capacity", 2])   )   


estimates.nonAfr_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Groups outside \n Sub-Saharan Africa", "Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa"),
  estimate = c(coef(nonAfr.lm.ucdp.2)["policing_subgroups"],
               coef(nonAfr.lm.ucdp.3)["policing_subgroups"],
               coef(nonAfr.lm.ucdp.4)["policing_subgroups"],
               coef(nonAfr.lm.ucdp.1)["policing_capacity"],
               coef(nonAfr.lm.ucdp.3)["policing_capacity"],
               coef(nonAfr.lm.ucdp.4)["policing_capacity"]),
  ci.lower = c(confint(nonAfr.lm.ucdp.2)["policing_subgroups", 1],
               confint(nonAfr.lm.ucdp.3)["policing_subgroups", 1],
               confint(nonAfr.lm.ucdp.4)["policing_subgroups", 1],
               confint(nonAfr.lm.ucdp.1)["policing_capacity", 1],
               confint(nonAfr.lm.ucdp.3)["policing_capacity", 1],
               confint(nonAfr.lm.ucdp.4)["policing_capacity", 1]),
  ci.upper = c(confint(nonAfr.lm.ucdp.2)["policing_subgroups", 2],
               confint(nonAfr.lm.ucdp.3)["policing_subgroups", 2],
               confint(nonAfr.lm.ucdp.4)["policing_subgroups", 2],
               confint(nonAfr.lm.ucdp.1)["policing_capacity", 2],
               confint(nonAfr.lm.ucdp.3)["policing_capacity", 2],
               confint(nonAfr.lm.ucdp.4)["policing_capacity", 2])   )  



## Interaction

# ACLED
model.ls <- list(
  nonAfr.inter.lmer.acled.1 <- lmer(formulae.me.inter.1[1, 1][[1]], 
                                    data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.inter.lmer.acled.2 <- lmer(formulae.me.inter.2[1, 1][[1]], 
                                    data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.inter.lmer.acled.3 <- lmer(formulae.me.inter.1[1, 2][[1]], 
                                    data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.inter.lmer.acled.4 <- lmer(formulae.me.inter.2[1, 2][[1]], 
                                    data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",])
)

# UCDP
model.ls <- list(
  nonAfr.inter.lmer.ucdp.1 <- lmer(formulae.me.inter.1[2, 1][[1]], 
                                   data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.inter.lmer.ucdp.2 <- lmer(formulae.me.inter.2[2, 1][[1]], 
                                   data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.inter.lmer.ucdp.3 <- lmer(formulae.me.inter.1[2, 2][[1]], 
                                   data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]),
  nonAfr.inter.lmer.ucdp.4 <- lmer(formulae.me.inter.2[2, 2][[1]], 
                                   data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",])
)



#Extract estimates
estimates.inter.nonAfr <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("Groups outside \n Sub-Saharan Africa", "Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa","Groups outside \n Sub-Saharan Africa"),
  estimate = c(coef(nonAfr.inter.lmer.acled.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(nonAfr.inter.lmer.acled.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(nonAfr.inter.lmer.acled.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(nonAfr.inter.lmer.acled.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(nonAfr.inter.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(nonAfr.inter.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(nonAfr.inter.lmer.ucdp.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(nonAfr.inter.lmer.ucdp.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(nonAfr.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(nonAfr.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(nonAfr.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(nonAfr.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(nonAfr.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(nonAfr.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(nonAfr.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(nonAfr.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(nonAfr.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(nonAfr.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(nonAfr.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(nonAfr.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(nonAfr.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(nonAfr.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(nonAfr.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(nonAfr.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )


estimates.nonAfr <- bind_rows(
  list(estimates.nonAfr_ucdp,
       estimates.nonAfr_acled,
       estimates.inter.nonAfr)
)


estimates.nonAfr$model_factor <- factor(estimates.nonAfr$model, 
                                        levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.nonAfr, file="estimates/estimates.nonAfr.RData" )


### > 2015-2020 ------

# ACLED
model.ls <- list(
  future.lm.acled.1 <- feols(formulae$twfe$twfe.1[[1]],
                             data = nc_tpi[nc_tpi$year > 2014,]),
  future.lm.acled.2 <- feols(formulae$twfe$twfe.2[[1]],
                             data = nc_tpi[nc_tpi$year > 2014,]),
  future.lm.acled.3 <- feols(formulae$twfe$twfe.3[[1]],
                             data = nc_tpi[nc_tpi$year > 2014,]),
  future.lm.acled.4 <- feols(formulae$twfe$twfe.4[[1]],
                             data = nc_tpi[nc_tpi$year > 2014,])
)


# UCDP
model.ls <- list(
  future.lm.ucdp.1 <- feols(formulae$twfe$twfe.1[[2]],
                            data = nc_tpi[nc_tpi$year > 2014,]),
  future.lm.ucdp.2 <- feols(formulae$twfe$twfe.2[[2]],
                            data = nc_tpi[nc_tpi$year > 2014,]),
  future.lm.ucdp.3 <- feols(formulae$twfe$twfe.3[[2]],
                            data = nc_tpi[nc_tpi$year > 2014,]),
  future.lm.ucdp.4 <- feols(formulae$twfe$twfe.4[[2]],
                            data = nc_tpi[nc_tpi$year > 2014,])
)



#Extract estimates

estimates.future_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("2015-2020", "2015-2020","2015-2020","2015-2020","2015-2020","2015-2020"),
  estimate = c(coef(future.lm.acled.2)["policing_subgroups"],
               coef(future.lm.acled.3)["policing_subgroups"],
               coef(future.lm.acled.4)["policing_subgroups"],
               coef(future.lm.acled.1)["policing_capacity"],
               coef(future.lm.acled.3)["policing_capacity"],
               coef(future.lm.acled.4)["policing_capacity"]),
  ci.lower = c(confint(future.lm.acled.2)["policing_subgroups", 1],
               confint(future.lm.acled.3)["policing_subgroups", 1],
               confint(future.lm.acled.4)["policing_subgroups", 1],
               confint(future.lm.acled.1)["policing_capacity", 1],
               confint(future.lm.acled.3)["policing_capacity", 1],
               confint(future.lm.acled.4)["policing_capacity", 1]),
  ci.upper = c(confint(future.lm.acled.2)["policing_subgroups", 2],
               confint(future.lm.acled.3)["policing_subgroups", 2],
               confint(future.lm.acled.4)["policing_subgroups", 2],
               confint(future.lm.acled.1)["policing_capacity", 2],
               confint(future.lm.acled.3)["policing_capacity", 2],
               confint(future.lm.acled.4)["policing_capacity", 2])   )   


estimates.future_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("2015-2020", "2015-2020","2015-2020","2015-2020","2015-2020","2015-2020"),
  estimate = c(coef(future.lm.ucdp.2)["policing_subgroups"],
               coef(future.lm.ucdp.3)["policing_subgroups"],
               coef(future.lm.ucdp.4)["policing_subgroups"],
               coef(future.lm.ucdp.1)["policing_capacity"],
               coef(future.lm.ucdp.3)["policing_capacity"],
               coef(future.lm.ucdp.4)["policing_capacity"]),
  ci.lower = c(confint(future.lm.ucdp.2)["policing_subgroups", 1],
               confint(future.lm.ucdp.3)["policing_subgroups", 1],
               confint(future.lm.ucdp.4)["policing_subgroups", 1],
               confint(future.lm.ucdp.1)["policing_capacity", 1],
               confint(future.lm.ucdp.3)["policing_capacity", 1],
               confint(future.lm.ucdp.4)["policing_capacity", 1]),
  ci.upper = c(confint(future.lm.ucdp.2)["policing_subgroups", 2],
               confint(future.lm.ucdp.3)["policing_subgroups", 2],
               confint(future.lm.ucdp.4)["policing_subgroups", 2],
               confint(future.lm.ucdp.1)["policing_capacity", 2],
               confint(future.lm.ucdp.3)["policing_capacity", 2],
               confint(future.lm.ucdp.4)["policing_capacity", 2])   )  


## Interaction

# ACLED
model.ls <- list(
  future.inter.lmer.acled.1 <- lmer(formulae.me.inter.1[1, 1][[1]], 
                                    data = nc_tpi[nc_tpi$year > 2014,]),
  future.inter.lmer.acled.2 <- lmer(formulae.me.inter.2[1, 1][[1]], 
                                    data = nc_tpi[nc_tpi$year > 2014,]),
  future.inter.lmer.acled.3 <- lmer(formulae.me.inter.1[1, 2][[1]], 
                                    data = nc_tpi[nc_tpi$year > 2014,]),
  future.inter.lmer.acled.4 <- lmer(formulae.me.inter.2[1, 2][[1]], 
                                    data = nc_tpi[nc_tpi$year > 2014,])
)


# UCDP
model.ls <- list(
  future.inter.lmer.ucdp.1 <- lmer(formulae.me.inter.1[2, 1][[1]], 
                                   data = nc_tpi[nc_tpi$year > 2014,]),
  future.inter.lmer.ucdp.2 <- lmer(formulae.me.inter.2[2, 1][[1]], 
                                   data = nc_tpi[nc_tpi$year > 2014,]),
  future.inter.lmer.ucdp.3 <- lmer(formulae.me.inter.1[2, 2][[1]], 
                                   data = nc_tpi[nc_tpi$year > 2014,]),
  future.inter.lmer.ucdp.4 <- lmer(formulae.me.inter.2[2, 2][[1]], 
                                   data = nc_tpi[nc_tpi$year > 2014,])
)


#Extract estimates
estimates.inter.future <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("2015-2020", "2015-2020","2015-2020","2015-2020","2015-2020","2015-2020","2015-2020","2015-2020"),
  estimate = c(coef(future.inter.lmer.acled.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(future.inter.lmer.acled.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(future.inter.lmer.acled.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(future.inter.lmer.acled.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(future.inter.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(future.inter.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(future.inter.lmer.ucdp.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(future.inter.lmer.ucdp.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(future.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(future.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(future.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(future.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(future.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(future.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(future.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(future.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(future.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(future.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(future.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(future.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(future.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(future.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(future.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(future.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )


estimates.future <- bind_rows(
  list(estimates.future_ucdp,
       estimates.future_acled,
       estimates.inter.future)
)



estimates.future$model_factor <- factor(estimates.future$model, 
                                        levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.future, file="estimates/estimates.future.RData" )




### > Dispute Resolution Only ----


#define new treatments
treat.dr <- c("disp_capacity", "policing_subgroups")
treat.dr.1 <- "disp_capacity"

cross.inter.dr.1 <- c("factor(two_class01_mcv)*disp_capacity", "policing_subgroups")


# Dispute Resolution
formulae$twfe$twfe.dr.1 <- lapply(outcomes_lags, construct_formula, treat = treat.dr.1, baseline = baseline.twfe)
formulae$twfe$twfe.dr.2 <- lapply(outcomes_lags, construct_formula, treat = treat.dr, baseline = baseline.twfe)
formulae$twfe$twfe.dr.3 <- lapply(outcomes_lags, construct_formula, treat = treat.dr, baseline = baseline.twfe, controls = controls.twfe)

formulae.me.dr.inter.1 <- sapply(cross.inter.dr.1, function(inter) {
  sapply(outcomes_lags, construct_formula_inter, inter = inter)
})

formulae.me.dr.inter.2 <- sapply(cross.inter.dr.1, function(inter) {
  sapply(outcomes_lags, construct_formula_inter_controls, inter = inter)
})



# ACLED
model.ls <- list(
  dr.lm.acled.1 <- feols(formulae$twfe$twfe.dr.1[[1]],
                         data = nc_tpi, vcov = "hetero"),
  dr.lm.acled.2 <- feols(formulae$twfe$twfe.dr.2[[1]],
                         data = nc_tpi, vcov = "hetero"),
  dr.lm.acled.3 <- feols(formulae$twfe$twfe.dr.3[[1]],
                         data = nc_tpi, vcov = "hetero")
)

# UCDP
model.ls <- list(
  dr.lm.ucdp.1 <- feols(formulae$twfe$twfe.dr.1[[2]],
                        data = nc_tpi, vcov = "hetero"),
  dr.lm.ucdp.2 <- feols(formulae$twfe$twfe.dr.2[[2]],
                        data = nc_tpi, vcov = "hetero"),
  dr.lm.ucdp.3 <- feols(formulae$twfe$twfe.dr.3[[2]],
                        data = nc_tpi, vcov = "hetero")
)




#Extract estimates

estimates.dr_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls"),
  variable = c("Supergroup Policing", "Supergroup Policing","Supergroup Policing"),
  specification = c("Dispute Resolution Only", "Dispute Resolution Only","Dispute Resolution Only"),
  estimate = c(coef(dr.lm.acled.1)["disp_capacity"],
               coef(dr.lm.acled.2)["disp_capacity"],
               coef(dr.lm.acled.3)["disp_capacity"]),
  ci.lower = c(confint(dr.lm.acled.1)["disp_capacity", 1],
               confint(dr.lm.acled.2)["disp_capacity", 1],
               confint(dr.lm.acled.3)["disp_capacity", 1]),
  ci.upper = c(confint(dr.lm.acled.1)["disp_capacity", 2],
               confint(dr.lm.acled.2)["disp_capacity", 2],
               confint(dr.lm.acled.3)["disp_capacity", 2])   )   


estimates.dr_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls"),
  variable = c("Supergroup Policing", "Supergroup Policing","Supergroup Policing"),
  specification = c("Dispute Resolution Only", "Dispute Resolution Only","Dispute Resolution Only"),
  estimate = c(coef(dr.lm.ucdp.1)["disp_capacity"],
               coef(dr.lm.ucdp.2)["disp_capacity"],
               coef(dr.lm.ucdp.3)["disp_capacity"]),
  ci.lower = c(confint(dr.lm.ucdp.1)["disp_capacity", 1],
               confint(dr.lm.ucdp.2)["disp_capacity", 1],
               confint(dr.lm.ucdp.3)["disp_capacity", 1]),
  ci.upper = c(confint(dr.lm.ucdp.1)["disp_capacity", 2],
               confint(dr.lm.ucdp.2)["disp_capacity", 2],
               confint(dr.lm.ucdp.3)["disp_capacity", 2])   )   




## Interaction

# ACLED
model.ls <- list(
  inter.dr.lmer.acled.1 <- lmer(formulae.me.dr.inter.1[1, 1][[1]], 
                                data = nc_tpi),
  inter.dr.lmer.acled.2 <- lmer(formulae.me.dr.inter.2[1, 1][[1]], 
                                data = nc_tpi)
)


# UCDP
model.ls <- list(
  inter.dr.lmer.ucdp.1 <- lmer(formulae.me.dr.inter.1[2, 1][[1]], 
                               data = nc_tpi),
  inter.dr.lmer.ucdp.2 <- lmer(formulae.me.dr.inter.2[2, 1][[1]], 
                               data = nc_tpi)
)



#Extract estimates
estimates.inter.dr <- data.frame(
  test = c("ACLED", "ACLED", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation"),
  specification = c("Dispute Resolution Only", "Dispute Resolution Only",
                    "Dispute Resolution Only","Dispute Resolution Only"),
  estimate = c(coef(inter.dr.lmer.acled.1)$cow[1,"factor(two_class01_mcv)1:disp_capacity"],
               coef(inter.dr.lmer.acled.2)$cow[1,"factor(two_class01_mcv)1:disp_capacity"],
               coef(inter.dr.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:disp_capacity"],
               coef(inter.dr.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:disp_capacity"]),
  ci.lower = c(confint(inter.dr.lmer.acled.1)["factor(two_class01_mcv)1:disp_capacity", 1],
               confint(inter.dr.lmer.acled.2)["factor(two_class01_mcv)1:disp_capacity", 1],
               confint(inter.dr.lmer.ucdp.1)["factor(two_class01_mcv)1:disp_capacity", 1],
               confint(inter.dr.lmer.ucdp.2)["factor(two_class01_mcv)1:disp_capacity", 1]),
  ci.upper = c(confint(inter.dr.lmer.acled.1)["factor(two_class01_mcv)1:disp_capacity", 2],
               confint(inter.dr.lmer.acled.2)["factor(two_class01_mcv)1:disp_capacity", 2],
               confint(inter.dr.lmer.ucdp.1)["factor(two_class01_mcv)1:disp_capacity", 2],
               confint(inter.dr.lmer.ucdp.2)["factor(two_class01_mcv)1:disp_capacity", 2]) )


estimates.dr <- bind_rows(
  list(estimates.dr_ucdp,
       estimates.dr_acled,
       estimates.inter.dr)
)


estimates.dr$model_factor <- factor(estimates.dr$model, 
                                    levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.dr, file="estimates/estimates.dr.RData" )



### > Security Only ----

treat.s <- c("sec_capacity", "policing_subgroups")
treat.s.1 <- "sec_capacity"
cross.inter.s.1 <- c("factor(two_class01_mcv)*sec_capacity", "policing_subgroups")

formulae$twfe$twfe.s.1 <- lapply(outcomes_lags, construct_formula, treat = treat.s.1, baseline = baseline.twfe)
formulae$twfe$twfe.s.2 <- lapply(outcomes_lags, construct_formula, treat = treat.s, baseline = baseline.twfe)
formulae$twfe$twfe.s.3 <- lapply(outcomes_lags, construct_formula, treat = treat.s, baseline = baseline.twfe, controls = controls.twfe)

formulae.me.s.inter.1 <- sapply(cross.inter.s.1, function(inter) {
  sapply(outcomes_lags, construct_formula_inter, inter = inter)
})

formulae.me.s.inter.2 <- sapply(cross.inter.s.1, function(inter) {
  sapply(outcomes_lags, construct_formula_inter_controls, inter = inter)
})



# ACLED
model.ls <- list(
  s.lm.acled.1 <- feols(formulae$twfe$twfe.s.1[[1]],
                        data = nc_tpi, vcov = "hetero"),
  s.lm.acled.2 <- feols(formulae$twfe$twfe.s.2[[1]],
                        data = nc_tpi, vcov = "hetero"),
  s.lm.acled.3 <- feols(formulae$twfe$twfe.s.3[[1]],
                        data = nc_tpi, vcov = "hetero")  
)


# UCDP
model.ls <- list(
  s.lm.ucdp.1 <- feols(formulae$twfe$twfe.s.1[[2]],
                       data = nc_tpi, vcov = "hetero"),
  s.lm.ucdp.2 <- feols(formulae$twfe$twfe.s.2[[2]],
                       data = nc_tpi, vcov = "hetero"),
  s.lm.ucdp.3 <- feols(formulae$twfe$twfe.s.3[[2]],
                       data = nc_tpi, vcov = "hetero")
)


#Extract estimates

estimates.s_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls"),
  variable = c("Supergroup Policing", "Supergroup Policing","Supergroup Policing"),
  specification = c("Security Only", "Security Only","Security Only"),
  estimate = c(coef(s.lm.acled.1)["sec_capacity"],
               coef(s.lm.acled.2)["sec_capacity"],
               coef(s.lm.acled.3)["sec_capacity"]),
  ci.lower = c(confint(s.lm.acled.1)["sec_capacity", 1],
               confint(s.lm.acled.2)["sec_capacity", 1],
               confint(s.lm.acled.3)["sec_capacity", 1]),
  ci.upper = c(confint(s.lm.acled.1)["sec_capacity", 2],
               confint(s.lm.acled.2)["sec_capacity", 2],
               confint(s.lm.acled.3)["sec_capacity", 2])   )   


estimates.s_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls"),
  variable = c("Supergroup Policing", "Supergroup Policing","Supergroup Policing"),
  specification = c("Security Only", "Security Only","Security Only"),
  estimate = c(coef(s.lm.ucdp.1)["sec_capacity"],
               coef(s.lm.ucdp.2)["sec_capacity"],
               coef(s.lm.ucdp.3)["sec_capacity"]),
  ci.lower = c(confint(s.lm.ucdp.1)["sec_capacity", 1],
               confint(s.lm.ucdp.2)["sec_capacity", 1],
               confint(s.lm.ucdp.3)["sec_capacity", 1]),
  ci.upper = c(confint(s.lm.ucdp.1)["sec_capacity", 2],
               confint(s.lm.ucdp.2)["sec_capacity", 2],
               confint(s.lm.ucdp.3)["sec_capacity", 2])   )   



## Interaction

# ACLED
model.ls <- list(
  inter.s.lmer.acled.1 <- lmer(formulae.me.s.inter.1[1, 1][[1]], 
                               data = nc_tpi),
  inter.s.lmer.acled.2 <- lmer(formulae.me.s.inter.2[1, 1][[1]], 
                               data = nc_tpi)
)

# UCDP
model.ls <- list(
  inter.s.lmer.ucdp.1 <- lmer(formulae.me.s.inter.1[2, 1][[1]], 
                              data = nc_tpi),
  inter.s.lmer.ucdp.2 <- lmer(formulae.me.s.inter.2[2, 1][[1]], 
                              data = nc_tpi)
)

#Extract estimates
estimates.inter.s <- data.frame(
  test = c("ACLED", "ACLED", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation"),
  specification = c("Security Only", "Security Only",
                    "Security Only","Security Only"),
  estimate = c(coef(inter.s.lmer.acled.1)$cow[1,"factor(two_class01_mcv)1:sec_capacity"],
               coef(inter.s.lmer.acled.2)$cow[1,"factor(two_class01_mcv)1:sec_capacity"],
               coef(inter.s.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:sec_capacity"],
               coef(inter.s.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:sec_capacity"]),
  ci.lower = c(confint(inter.s.lmer.acled.1)["factor(two_class01_mcv)1:sec_capacity", 1],
               confint(inter.s.lmer.acled.2)["factor(two_class01_mcv)1:sec_capacity", 1],
               confint(inter.s.lmer.ucdp.1)["factor(two_class01_mcv)1:sec_capacity", 1],
               confint(inter.s.lmer.ucdp.2)["factor(two_class01_mcv)1:sec_capacity", 1]),
  ci.upper = c(confint(inter.s.lmer.acled.1)["factor(two_class01_mcv)1:sec_capacity", 2],
               confint(inter.s.lmer.acled.2)["factor(two_class01_mcv)1:sec_capacity", 2],
               confint(inter.s.lmer.ucdp.1)["factor(two_class01_mcv)1:sec_capacity", 2],
               confint(inter.s.lmer.ucdp.2)["factor(two_class01_mcv)1:sec_capacity", 2]) )


estimates.s <- bind_rows(
  list(estimates.s_ucdp,
       estimates.s_acled,
       estimates.inter.s)
)

estimates.s$model_factor <- factor(estimates.s$model, 
                                   levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.s, file="estimates/estimates.s.RData" )



### > Controlling for All Subgroups ----

# ACLED
model.ls <- list(
  amar.lm.acled.1 <- feols(formulae$twfe$twfe.1[[7]],
                           data = nc_tpi, vcov = "hetero"),
  amar.lm.acled.2 <- feols(formulae$twfe$twfe.2[[7]],
                           data = nc_tpi, vcov = "hetero"),
  amar.lm.acled.3 <- feols(formulae$twfe$twfe.3[[7]],
                           data = nc_tpi, vcov = "hetero"),
  amar.lm.acled.4 <- feols(formulae$twfe$twfe.4[[7]],
                           data = nc_tpi, vcov = "hetero")
)


#Extract estimates

estimates.subct_acled <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "ACLED", "ACLED", "ACLED"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Controlling for All Supgroups", "Controlling for All Supgroups","Controlling for All Supgroups","Controlling for All Supgroups","Controlling for All Supgroups","Controlling for All Supgroups"),
  estimate = c(coef(amar.lm.acled.2)["policing_subgroups"],
               coef(amar.lm.acled.3)["policing_subgroups"],
               coef(amar.lm.acled.4)["policing_subgroups"],
               coef(amar.lm.acled.1)["policing_capacity"],
               coef(amar.lm.acled.3)["policing_capacity"],
               coef(amar.lm.acled.4)["policing_capacity"]),
  ci.lower = c(confint(amar.lm.acled.2)["policing_subgroups", 1],
               confint(amar.lm.acled.3)["policing_subgroups", 1],
               confint(amar.lm.acled.4)["policing_subgroups", 1],
               confint(amar.lm.acled.1)["policing_capacity", 1],
               confint(amar.lm.acled.3)["policing_capacity", 1],
               confint(amar.lm.acled.4)["policing_capacity", 1]),
  ci.upper = c(confint(amar.lm.acled.2)["policing_subgroups", 2],
               confint(amar.lm.acled.3)["policing_subgroups", 2],
               confint(amar.lm.acled.4)["policing_subgroups", 2],
               confint(amar.lm.acled.1)["policing_capacity", 2],
               confint(amar.lm.acled.3)["policing_capacity", 2],
               confint(amar.lm.acled.4)["policing_capacity", 2])   )   


# UCDP
model.ls <- list(
  amar.lm.ucdp.1 <- feols(formulae$twfe$twfe.1[[8]],
                          data = nc_tpi, vcov = "hetero"),
  amar.lm.ucdp.2 <- feols(formulae$twfe$twfe.2[[8]],
                          data = nc_tpi, vcov = "hetero"),
  amar.lm.ucdp.3 <- feols(formulae$twfe$twfe.3[[8]],
                          data = nc_tpi, vcov = "hetero"),
  amar.lm.ucdp.4 <- feols(formulae$twfe$twfe.4[[8]],
                          data = nc_tpi, vcov = "hetero")
)


estimates.subct_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", 
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "Supergroup Policing",
               "Supergroup Policing","Supergroup Policing"),
  specification = c("Controlling for All Supgroups", "Controlling for All Supgroups","Controlling for All Supgroups","Controlling for All Supgroups","Controlling for All Supgroups","Controlling for All Supgroups"),
  estimate = c(coef(amar.lm.ucdp.2)["policing_subgroups"],
               coef(amar.lm.ucdp.3)["policing_subgroups"],
               coef(amar.lm.ucdp.4)["policing_subgroups"],
               coef(amar.lm.ucdp.1)["policing_capacity"],
               coef(amar.lm.ucdp.3)["policing_capacity"],
               coef(amar.lm.ucdp.4)["policing_capacity"]),
  ci.lower = c(confint(amar.lm.ucdp.2)["policing_subgroups", 1],
               confint(amar.lm.ucdp.3)["policing_subgroups", 1],
               confint(amar.lm.ucdp.4)["policing_subgroups", 1],
               confint(amar.lm.ucdp.1)["policing_subgroups", 1],
               confint(amar.lm.ucdp.3)["policing_capacity", 1],
               confint(amar.lm.ucdp.4)["policing_capacity", 1]),
  ci.upper = c(confint(amar.lm.ucdp.2)["policing_subgroups", 2],
               confint(amar.lm.ucdp.3)["policing_subgroups", 2],
               confint(amar.lm.ucdp.4)["policing_subgroups", 2],
               confint(amar.lm.ucdp.1)["policing_capacity", 2],
               confint(amar.lm.ucdp.3)["policing_capacity", 2],
               confint(amar.lm.ucdp.4)["policing_capacity", 2])   )  



## Interaction
# ACLED
model.ls <- list(
  amar.inter.lmer.acled.1 <- lmer(formulae.me.inter.1[7, 1][[1]], 
                                  data = nc_tpi),
  amar.inter.lmer.acled.2 <- lmer(formulae.me.inter.2[7, 1][[1]], 
                                  data = nc_tpi),
  amar.inter.lmer.acled.3 <- lmer(formulae.me.inter.1[7, 2][[1]], 
                                  data = nc_tpi),
  amar.inter.lmer.acled.4 <- lmer(formulae.me.inter.2[7, 2][[1]], 
                                  data = nc_tpi)
)


# UCDP
model.ls <- list(
  amar.inter.lmer.ucdp.1 <- lmer(formulae.me.inter.1[8, 1][[1]], 
                                 data = nc_tpi),
  amar.inter.lmer.ucdp.2 <- lmer(formulae.me.inter.2[8, 1][[1]], 
                                 data = nc_tpi),
  amar.inter.lmer.ucdp.3 <- lmer(formulae.me.inter.1[8, 2][[1]], 
                                 data = nc_tpi),
  amar.inter.lmer.ucdp.4 <- lmer(formulae.me.inter.2[8, 2][[1]], 
                                 data = nc_tpi)
)



#Extract estimates
estimates.inter.subct <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("Controlling for All Supgroups", "Controlling for All Supgroups","Controlling for All Supgroups",
                    "Controlling for All Supgroups","Controlling for All Supgroups","Controlling for All Supgroups",
                    "Controlling for All Supgroups","Controlling for All Supgroups"),
  estimate = c(coef(amar.inter.lmer.acled.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(amar.inter.lmer.acled.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(amar.inter.lmer.acled.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(amar.inter.lmer.acled.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(amar.inter.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(amar.inter.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(amar.inter.lmer.ucdp.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(amar.inter.lmer.ucdp.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(amar.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(amar.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(amar.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(amar.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(amar.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(amar.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(amar.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(amar.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(amar.inter.lmer.acled.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(amar.inter.lmer.acled.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(amar.inter.lmer.acled.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(amar.inter.lmer.acled.4)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(amar.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(amar.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(amar.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(amar.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )


estimates.subct <- bind_rows(
  list(estimates.subct_ucdp,
       estimates.subct_acled,
       estimates.inter.subct)
)


estimates.subct$model_factor <- factor(estimates.subct$model, 
                                       levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.subct, file="estimates/estimates.subct.RData" )


### > Controlling for Jurisdictional Conflict ----

# See above: Mediation I



### > Territory Conflicts Only -----

model.ls <- list(
  terr.lm.1 <- feols(formulae$twfe$twfe.1[[3]],
                  data = nc_tpi, vcov = "hetero"),
  terr.lm.2 <- feols(formulae$twfe$twfe.2[[3]],
                  data = nc_tpi, vcov = "hetero"),
  terr.lm.3 <- feols(formulae$twfe$twfe.3[[3]],
                  data = nc_tpi, vcov = "hetero"),
  terr.lm.4 <- feols(formulae$twfe$twfe.4[[3]],
                  data = nc_tpi, vcov = "hetero")
)



## Interaction
model.ls <- list(
  terr.inter.lmer.1 <- lmer(formulae.me.inter.1[3, 1][[1]], 
                            data = nc_tpi),
  terr.inter.lmer.2 <- lmer(formulae.me.inter.2[3, 1][[1]], 
                            data = nc_tpi),
  terr.inter.lmer.3 <- lmer(formulae.me.inter.1[3, 2][[1]], 
                            data = nc_tpi),
  terr.inter.lmer.4 <- lmer(formulae.me.inter.2[3, 2][[1]], 
                            data = nc_tpi)
)


#Extract estimates
estimates.terr_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", 
               "Supergroup Policing", "Supergroup Policing","Supergroup Policing"),
  specification = c("Territory Conflicts Only", "Territory Conflicts Only", "Territory Conflicts Only",
                    "Territory Conflicts Only","Territory Conflicts Only","Territory Conflicts Only"),
  estimate = c(coef(terr.lm.2)["policing_subgroups"],
               coef(terr.lm.3)["policing_subgroups"],
               coef(terr.lm.4)["policing_subgroups"],
               coef(terr.lm.1)["policing_capacity"],
               coef(terr.lm.3)["policing_capacity"],
               coef(terr.lm.4)["policing_capacity"]),
  ci.lower = c(confint(terr.lm.2)["policing_subgroups", 1],
               confint(terr.lm.3)["policing_subgroups", 1],
               confint(terr.lm.4)["policing_subgroups", 1],
               confint(terr.lm.1)["policing_capacity", 1],
               confint(terr.lm.3)["policing_capacity", 1],
               confint(terr.lm.4)["policing_capacity", 1]),
  ci.upper = c(confint(terr.lm.2)["policing_subgroups", 2],
               confint(terr.lm.3)["policing_subgroups", 2],
               confint(terr.lm.4)["policing_subgroups", 2],
               confint(terr.lm.1)["policing_capacity", 2],
               confint(terr.lm.3)["policing_capacity", 2],
               confint(terr.lm.4)["policing_capacity", 2])   )  


estimates.inter.terr <- data.frame(
  test = c("UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("Territory Conflicts Only","Territory Conflicts Only","Territory Conflicts Only","Territory Conflicts Only"),
  estimate = c(coef(terr.inter.lmer.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(terr.inter.lmer.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(terr.inter.lmer.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(terr.inter.lmer.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(terr.inter.lmer.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(terr.inter.lmer.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(terr.inter.lmer.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(terr.inter.lmer.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(terr.inter.lmer.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(terr.inter.lmer.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(terr.inter.lmer.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(terr.inter.lmer.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )


estimates.terr <- bind_rows(
  list(estimates.terr_ucdp,
       estimates.inter.terr)
)

estimates.terr$model_factor <- factor(estimates.terr$model, 
                                 levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.terr, file="estimates/estimates.terr.RData" )



### > Authority Conflicts Only -----

model.ls <- list(
  aut.lm.1 <- feols(formulae$twfe$twfe.1[[4]],
                 data = nc_tpi, vcov = "hetero"),
  aut.lm.2 <- feols(formulae$twfe$twfe.2[[4]],
                 data = nc_tpi, vcov = "hetero"),
  aut.lm.3 <- feols(formulae$twfe$twfe.3[[4]],
                 data = nc_tpi, vcov = "hetero"),
  aut.lm.4 <- feols(formulae$twfe$twfe.4[[4]],
                 data = nc_tpi, vcov = "hetero")
)


## Interaction
model.ls <- list(
  aut.inter.lmer.1 <- lmer(formulae.me.inter.1[4, 1][[1]], 
                           data = nc_tpi),
  aut.inter.lmer.2 <- lmer(formulae.me.inter.2[4, 1][[1]], 
                           data = nc_tpi),
  aut.inter.lmer.3 <- lmer(formulae.me.inter.1[4, 2][[1]], 
                           data = nc_tpi),
  aut.inter.lmer.4 <- lmer(formulae.me.inter.2[4, 2][[1]], 
                           data = nc_tpi)
)

#Extract estimates
estimates.aut_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions",  "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", 
               "Supergroup Policing", "Supergroup Policing","Supergroup Policing"),
  specification = c("Authority Conflicts Only", "Authority Conflicts Only","Authority Conflicts Only",
                    "Authority Conflicts Only","Authority Conflicts Only","Authority Conflicts Only"),
  estimate = c(coef(aut.lm.2)["policing_subgroups"],
               coef(aut.lm.3)["policing_subgroups"],
               coef(aut.lm.4)["policing_subgroups"],
               coef(aut.lm.1)["policing_capacity"],
               coef(aut.lm.3)["policing_capacity"],
               coef(aut.lm.4)["policing_capacity"]),
  ci.lower = c(confint(aut.lm.2)["policing_subgroups", 1],
               confint(aut.lm.3)["policing_subgroups", 1],
               confint(aut.lm.4)["policing_subgroups", 1],
               confint(aut.lm.1)["policing_capacity", 1],
               confint(aut.lm.3)["policing_capacity", 1],
               confint(aut.lm.4)["policing_capacity", 1]),
  ci.upper = c(confint(aut.lm.2)["policing_subgroups", 2],
               confint(aut.lm.3)["policing_subgroups", 2],
               confint(aut.lm.4)["policing_subgroups", 2],
               confint(aut.lm.1)["policing_capacity", 2],
               confint(aut.lm.3)["policing_capacity", 2],
               confint(aut.lm.4)["policing_capacity", 2])   )  


estimates.inter.aut <- data.frame(
  test = c("UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("Authority Conflicts Only","Authority Conflicts Only","Authority Conflicts Only","Authority Conflicts Only"),
  estimate = c(coef(aut.inter.lmer.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(aut.inter.lmer.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(aut.inter.lmer.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(aut.inter.lmer.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(aut.inter.lmer.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(aut.inter.lmer.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(aut.inter.lmer.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(aut.inter.lmer.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(aut.inter.lmer.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(aut.inter.lmer.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(aut.inter.lmer.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(aut.inter.lmer.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )


estimates.aut <- bind_rows(
  list(estimates.aut_ucdp,
       estimates.inter.aut)
)

estimates.aut$model_factor <- factor(estimates.aut$model, 
                                      levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.aut, file="estimates/estimates.aut.RData" )



### > Non-Geo Outcome ----

# UCDP
model.ls <- list(
  ng.lm.ucdp.1 <- feols(formulae$twfe$twfe.1[[9]],
                    data = nc_tpi, vcov = "hetero"),
  ng.lm.ucdp.2 <- feols(formulae$twfe$twfe.2[[9]],
                    data = nc_tpi, vcov = "hetero"),
  ng.lm.ucdp.3 <- feols(formulae$twfe$twfe.3[[9]],
                    data = nc_tpi, vcov = "hetero"),
  ng.lm.ucdp.4 <- feols(formulae$twfe$twfe.4[[9]],
                    data = nc_tpi, vcov = "hetero")
) 

estimates.ng_ucdp <- data.frame(
  test = c( "UCDP", "UCDP", "UCDP",
            "UCDP", "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls", "One treatment", "Both treatments", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions", 
               "Supergroup Policing","Supergroup Policing","Supergroup Policing"),
  specification = c("Non-Geo Outcome", "Non-Geo Outcome","Non-Geo Outcome",
                    "Non-Geo Outcome", "Non-Geo Outcome","Non-Geo Outcome"),
  estimate = c(coef(ng.lm.ucdp.2)["policing_subgroups"],
               coef(ng.lm.ucdp.3)["policing_subgroups"],
               coef(ng.lm.ucdp.4)["policing_subgroups"],
               coef(ng.lm.ucdp.1)["policing_capacity"],
               coef(ng.lm.ucdp.3)["policing_capacity"],
               coef(ng.lm.ucdp.4)["policing_capacity"]),
  ci.lower = c(confint(ng.lm.ucdp.2)["policing_subgroups", 1],
               confint(ng.lm.ucdp.3)["policing_subgroups", 1],
               confint(ng.lm.ucdp.4)["policing_subgroups", 1],
               confint(ng.lm.ucdp.1)["policing_capacity", 1],
               confint(ng.lm.ucdp.3)["policing_capacity", 1],
               confint(ng.lm.ucdp.4)["policing_capacity", 1]),
  ci.upper = c(confint(ng.lm.ucdp.2)["policing_subgroups", 2],
               confint(ng.lm.ucdp.3)["policing_subgroups", 2],
               confint(ng.lm.ucdp.4)["policing_subgroups", 2],
               confint(ng.lm.ucdp.1)["policing_capacity", 2],
               confint(ng.lm.ucdp.3)["policing_capacity", 2],
               confint(ng.lm.ucdp.4)["policing_capacity", 2])   )  


# Interaction
model.ls <- list(
  ng.inter.lmer.ucdp.1 <- lmer(formulae.me.inter.1[9, 1][[1]], 
                           data = nc_tpi),
  ng.inter.lmer.ucdp.2 <- lmer(formulae.me.inter.2[9, 1][[1]], 
                           data = nc_tpi),
  ng.inter.lmer.ucdp.3 <- lmer(formulae.me.inter.1[9, 2][[1]], 
                           data = nc_tpi),
  ng.inter.lmer.ucdp.4 <- lmer(formulae.me.inter.2[9, 2][[1]], 
                           data = nc_tpi)
)



#Extract estimates
estimates.inter.ng <- data.frame(
  test = c("UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation", 
               "Supergroup Policing x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation"),
  specification = c("Non-Geo Outcome","Non-Geo Outcome","Non-Geo Outcome","Non-Geo Outcome"),
  estimate = c(coef(ng.inter.lmer.ucdp.1)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(ng.inter.lmer.ucdp.2)$cow[1,"factor(two_class01_mcv)1:policing_capacity"],
               coef(ng.inter.lmer.ucdp.3)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"],
               coef(ng.inter.lmer.ucdp.4)$cow[1,"factor(two_class01_mcv)1:policing_subgroups"]),
  ci.lower = c(confint(ng.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(ng.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 1],
               confint(ng.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 1],
               confint(ng.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 1]),
  ci.upper = c(confint(ng.inter.lmer.ucdp.1)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(ng.inter.lmer.ucdp.2)["factor(two_class01_mcv)1:policing_capacity", 2],
               confint(ng.inter.lmer.ucdp.3)["factor(two_class01_mcv)1:policing_subgroups", 2],
               confint(ng.inter.lmer.ucdp.4)["factor(two_class01_mcv)1:policing_subgroups", 2]) )



estimates.ng <- bind_rows(
  list(estimates.ng_ucdp,
       estimates.inter.ng)
)

estimates.ng$model_factor <- factor(estimates.ng$model, 
                                     levels = c("Incl. controls", "Both treatments", "One treatment"))
save(estimates.ng, file="estimates/estimates.ng.RData" )



### > Constitutional Regulation: Min-Max ------

# ACLED
model.ls <- list(
  inter.lmer.rega.acled.1 <- lmer(formulae.me.inter.1[1, 7][[1]], 
                                  data = nc_tpi),
  inter.lmer.rega.acled.2 <- lmer(formulae.me.inter.2[1, 7][[1]], 
                                  data = nc_tpi),
  inter.lmer.rega.acled.3 <- lmer(formulae.me.inter.1[1, 8][[1]], 
                                  data = nc_tpi),
  inter.lmer.rega.acled.4 <- lmer(formulae.me.inter.2[1, 8][[1]], 
                                  data = nc_tpi)
)

# UCDP
model.ls <- list(
  inter.lmer.rega.ucdp.1 <- lmer(formulae.me.inter.1[2, 7][[1]], 
                                 data = nc_tpi),
  inter.lmer.rega.ucdp.2 <- lmer(formulae.me.inter.2[2, 7][[1]], 
                                 data = nc_tpi),
  inter.lmer.rega.ucdp.3 <- lmer(formulae.me.inter.1[2, 8][[1]], 
                                 data = nc_tpi),
  inter.lmer.rega.ucdp.4 <- lmer(formulae.me.inter.2[2, 8][[1]], 
                                 data = nc_tpi)
)



#Extract estimates
estimates.inter.rega <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation (min-max)", 
               "Supergroup Policing x \n Constitutional Regulation (min-max)",
               "No. of Subgroups x \n Constitutional Regulation (min-max)",
               "No. of Subgroups x \n Constitutional Regulation (min-max)"),
  specification = c("Constitutional Regulation: Min-Max", "Constitutional Regulation: Min-Max","Constitutional Regulation: Min-Max","Constitutional Regulation: Min-Max","Constitutional Regulation: Min-Max","Constitutional Regulation: Min-Max","Constitutional Regulation: Min-Max","Constitutional Regulation: Min-Max"),
  estimate = c(coef(inter.lmer.rega.acled.1)$cow[1,"REGA_n_mcv:policing_capacity"],
               coef(inter.lmer.rega.acled.2)$cow[1,"REGA_n_mcv:policing_capacity"],
               coef(inter.lmer.rega.acled.3)$cow[1,"REGA_n_mcv:policing_subgroups"],
               coef(inter.lmer.rega.acled.4)$cow[1,"REGA_n_mcv:policing_subgroups"],
               coef(inter.lmer.rega.ucdp.1)$cow[1,"REGA_n_mcv:policing_capacity"],
               coef(inter.lmer.rega.ucdp.2)$cow[1,"REGA_n_mcv:policing_capacity"],
               coef(inter.lmer.rega.ucdp.3)$cow[1,"REGA_n_mcv:policing_subgroups"],
               coef(inter.lmer.rega.ucdp.4)$cow[1,"REGA_n_mcv:policing_subgroups"]),
  ci.lower = c(confint(inter.lmer.rega.acled.1)["REGA_n_mcv:policing_capacity", 1],
               confint(inter.lmer.rega.acled.2)["REGA_n_mcv:policing_capacity", 1],
               confint(inter.lmer.rega.acled.3)["REGA_n_mcv:policing_subgroups", 1],
               confint(inter.lmer.rega.acled.4)["REGA_n_mcv:policing_subgroups", 1],
               confint(inter.lmer.rega.ucdp.1)["REGA_n_mcv:policing_capacity", 1],
               confint(inter.lmer.rega.ucdp.2)["REGA_n_mcv:policing_capacity", 1],
               confint(inter.lmer.rega.ucdp.3)["REGA_n_mcv:policing_subgroups", 1],
               confint(inter.lmer.rega.ucdp.4)["REGA_n_mcv:policing_subgroups", 1]),
  ci.upper = c(confint(inter.lmer.rega.acled.1)["REGA_n_mcv:policing_capacity", 2],
               confint(inter.lmer.rega.acled.2)["REGA_n_mcv:policing_capacity", 2],
               confint(inter.lmer.rega.acled.3)["REGA_n_mcv:policing_subgroups", 2],
               confint(inter.lmer.rega.acled.4)["REGA_n_mcv:policing_subgroups", 2],
               confint(inter.lmer.rega.ucdp.1)["REGA_n_mcv:policing_capacity", 2],
               confint(inter.lmer.rega.ucdp.2)["REGA_n_mcv:policing_capacity", 2],
               confint(inter.lmer.rega.ucdp.3)["REGA_n_mcv:policing_subgroups", 2],
               confint(inter.lmer.rega.ucdp.4)["REGA_n_mcv:policing_subgroups", 2]) )


estimates.inter.rega$model_factor <- factor(estimates.inter.rega$model, 
                                            levels = c("Incl. controls", "Excl. controls"))


save(estimates.inter.rega, file="estimates/estimates.inter.rega.RData" )


# >> Figure A5 ----

interplot(m=inter.lmer.rega.acled.2, var1="policing_capacity", var2="REGA_n_mcv",  hist=T,   ci=0.95) +
  ylab("Marginal Effect of Supergroup on Communal Conflict (ACLED)") +
  xlab("Constitutional Regulation (min-max)") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=16)) +
  theme(axis.title.y = element_text(size=16)) +
  theme(plot.title = element_text(face="bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  aes(color = "pink") + theme(legend.position="none") 


ggsave("interplot_cap_acled.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =18,
       units = "cm",
       dpi=700)


interplot(m=inter.lmer.rega.ucdp.2, var1="policing_capacity", var2="REGA_n_mcv",  hist=T,   ci=0.95) +
  ylab("Marginal Effect of Supergroup on Communal Conflict (UCDP)") +
  xlab("Constitutional Regulation (min-max)") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=16)) +
  theme(axis.title.y = element_text(size=16)) +
  theme(plot.title = element_text(face="bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  aes(color = "pink") + theme(legend.position="none") 


ggsave("interplot_cap_ucdp.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =18,
       units = "cm",
       dpi=700)


interplot(m=inter.lmer.rega.acled.4, var1="policing_subgroups", var2="REGA_n_mcv",  hist=T,   ci=0.95,) +
  ylab("Marginal Effect of Subgroups on Communal Conflict (ACLED)") +
  xlab("Constitutional Regulation (min-max)") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=16)) +
  theme(axis.title.y = element_text(size=16)) +
  theme(plot.title = element_text(face="bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  aes(color = "pink") + theme(legend.position="none") 


ggsave("interplot_horiz_acled.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =18,
       units = "cm",
       dpi=700)

interplot(m=inter.lmer.rega.ucdp.4, var1="policing_subgroups", var2="REGA_n_mcv",  hist=T,   ci=0.95) +
  ylab("Marginal Effect of Subgroups on Communal Conflict (UCDP)") +
  xlab("Constitutional Regulation (min-max)") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=16)) +
  theme(axis.title.y = element_text(size=16)) +
  theme(plot.title = element_text(face="bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  aes(color = "pink") + theme(legend.position="none") 


ggsave("interplot_horiz_ucdp.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =18,
       units = "cm",
       dpi=700)




### > Constitutional Regulation: PC1 ----

# ACLED
model.ls <- list(
  inter.lmer.pc1.acled.1 <- lmer(formulae.me.inter.1[1, 5][[1]], 
                           data = nc_tpi),
  inter.lmer.pc1.acled.2 <- lmer(formulae.me.inter.2[1, 5][[1]], 
                           data = nc_tpi),
  inter.lmer.pc1.acled.3 <- lmer(formulae.me.inter.1[1, 6][[1]], 
                           data = nc_tpi),
  inter.lmer.pc1.acled.4 <- lmer(formulae.me.inter.2[1, 6][[1]], 
                           data = nc_tpi)
)


# UCDP
model.ls <- list(
  inter.lmer.pc1.ucdp.1 <- lmer(formulae.me.inter.1[2, 5][[1]], 
                          data = nc_tpi),
  inter.lmer.pc1.ucdp.2 <- lmer(formulae.me.inter.2[2, 5][[1]], 
                          data = nc_tpi),
  inter.lmer.pc1.ucdp.3 <- lmer(formulae.me.inter.1[2, 6][[1]], 
                          data = nc_tpi),
  inter.lmer.pc1.ucdp.4 <- lmer(formulae.me.inter.2[2, 6][[1]], 
                          data = nc_tpi)
)



#Extract estimates
estimates.inter.pc1 <- data.frame(
  test = c("ACLED", "ACLED", "ACLED", 
           "ACLED", "UCDP", "UCDP", 
           "UCDP", "UCDP"),
  model = c("Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("Supergroup Policing x \n Constitutional Regulation (First Principal Component)", 
               "Supergroup Policing x \n Constitutional Regulation (First Principal Component)",
               "No. of Subgroups x \n Constitutional Regulation (First Principal Component)",
               "No. of Subgroups x \n Constitutional Regulation (First Principal Component)"),
  specification = c("Constitutional Regulation: PC1", "Constitutional Regulation: PC1","Constitutional Regulation: PC1","Constitutional Regulation: PC1","Constitutional Regulation: PC1","Constitutional Regulation: PC1","Constitutional Regulation: PC1","Constitutional Regulation: PC1"),
  estimate = c(coef(inter.lmer.pc1.acled.1)$cow[1,"pc1_mcv:policing_capacity"],
               coef(inter.lmer.pc1.acled.2)$cow[1,"pc1_mcv:policing_capacity"],
               coef(inter.lmer.pc1.acled.3)$cow[1,"pc1_mcv:policing_subgroups"],
               coef(inter.lmer.pc1.acled.4)$cow[1,"pc1_mcv:policing_subgroups"],
               coef(inter.lmer.pc1.ucdp.1)$cow[1,"pc1_mcv:policing_capacity"],
               coef(inter.lmer.pc1.ucdp.2)$cow[1,"pc1_mcv:policing_capacity"],
               coef(inter.lmer.pc1.ucdp.3)$cow[1,"pc1_mcv:policing_subgroups"],
               coef(inter.lmer.pc1.ucdp.4)$cow[1,"pc1_mcv:policing_subgroups"]),
  ci.lower = c(confint(inter.lmer.pc1.acled.1)["pc1_mcv:policing_capacity", 1],
               confint(inter.lmer.pc1.acled.2)["pc1_mcv:policing_capacity", 1],
               confint(inter.lmer.pc1.acled.3)["pc1_mcv:policing_subgroups", 1],
               confint(inter.lmer.pc1.acled.4)["pc1_mcv:policing_subgroups", 1],
               confint(inter.lmer.pc1.ucdp.1)["pc1_mcv:policing_capacity", 1],
               confint(inter.lmer.pc1.ucdp.2)["pc1_mcv:policing_capacity", 1],
               confint(inter.lmer.pc1.ucdp.3)["pc1_mcv:policing_subgroups", 1],
               confint(inter.lmer.pc1.ucdp.4)["pc1_mcv:policing_subgroups", 1]),
  ci.upper = c(confint(inter.lmer.pc1.acled.1)["pc1_mcv:policing_capacity", 2],
               confint(inter.lmer.pc1.acled.2)["pc1_mcv:policing_capacity", 2],
               confint(inter.lmer.pc1.acled.3)["pc1_mcv:policing_subgroups", 2],
               confint(inter.lmer.pc1.acled.4)["pc1_mcv:policing_subgroups", 2],
               confint(inter.lmer.pc1.ucdp.1)["pc1_mcv:policing_capacity", 2],
               confint(inter.lmer.pc1.ucdp.2)["pc1_mcv:policing_capacity", 2],
               confint(inter.lmer.pc1.ucdp.3)["pc1_mcv:policing_subgroups", 2],
               confint(inter.lmer.pc1.ucdp.4)["pc1_mcv:policing_subgroups", 2]) )


estimates.inter.pc1$model_factor <- factor(estimates.inter.pc1$model, 
                                           levels = c("Incl. controls", "Excl. controls"))

save(estimates.inter.pc1, file="estimates/estimates.inter.pc1.RData" )


# >> Figure A4 ----

interplot(m=inter.lmer.pc1.acled.2, var1="policing_capacity", var2="pc1_mcv",  hist=T,   ci=0.95) +
  ylab("Marginal Effect of Supergroup on Communal Conflict (ACLED)") +
  xlab("Constitutional Regulation (PC1)") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=16)) +
  theme(axis.title.y = element_text(size=16)) +
  theme(plot.title = element_text(face="bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  aes(color = "pink") + theme(legend.position="none") 


ggsave("interplot_pc1_cap_acled.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =18,
       units = "cm",
       dpi=700)


interplot(m=inter.lmer.pc1.ucdp.2, var1="policing_capacity", var2="pc1_mcv",  hist=T,   ci=0.95) +
  ylab("Marginal Effect of Supergroup on Communal Conflict (UCDP)") +
  xlab("Constitutional Regulation (PC1)") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=16)) +
  theme(axis.title.y = element_text(size=16)) +
  theme(plot.title = element_text(face="bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  aes(color = "pink") + theme(legend.position="none") 


ggsave("interplot_pc1_cap_ucdp.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =18,
       units = "cm",
       dpi=700)


interplot(m=inter.lmer.pc1.acled.4, var1="policing_subgroups", var2="pc1_mcv",  hist=T,   ci=0.95,) +
  ylab("Marginal Effect of Subgroups on Communal Conflict (ACLED)") +
  xlab("Constitutional Regulation (PC1)") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=16)) +
  theme(axis.title.y = element_text(size=16)) +
  theme(plot.title = element_text(face="bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  aes(color = "pink") + theme(legend.position="none") 


ggsave("interplot_pc1_horiz_acled.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =18,
       units = "cm",
       dpi=700)

interplot(m=inter.lmer.pc1.ucdp.4, var1="policing_subgroups", var2="pc1_mcv",  hist=T,   ci=0.95) +
  ylab("Marginal Effect of Subgroups on Communal Conflict (UCDP)") +
  xlab("Constitutional Regulation (PC1)") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=16)) +
  theme(axis.title.y = element_text(size=16)) +
  theme(plot.title = element_text(face="bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  aes(color = "pink") + theme(legend.position="none") 


ggsave("interplot_pc1_horiz_ucdp.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =18,
       units = "cm",
       dpi=700)






### > Plotting Figure A3 ---- 


#see specifications above

estimates.robust <- bind_rows(
  list(
    estimates,
    estimates.inter,
    estimates.lmer,
    estimates.future,
    estimates.all_grps,
    estimates.nonAfr,
    estimates.dr,
    estimates.s,
    estimates.med,
    estimates.inter.pc1,
    estimates.inter.rega,
    estimates.aut,
    estimates.terr,
    estimates.ng,
    estimates.subct
    
  )
)

save(estimates.robust, file= "estimates/estimates.robustness.RData")


#adding NAs for models without associated test for side by side presentation in graph 
estimates.dr.na <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "UCDP", "UCDP", "UCDP",
            "ACLED", "ACLED", 
            "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls",
            "One treatment", "Both treatments", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions",
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions",
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions",
               "No. of Subgroups x \n Constitutional Regulation","No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation","No. of Subgroups x \n Constitutional Regulation" ),
  specification = c("Dispute Resolution Only", "Dispute Resolution Only","Dispute Resolution Only",
                    "Dispute Resolution Only", "Dispute Resolution Only","Dispute Resolution Only",
                    "Dispute Resolution Only", "Dispute Resolution Only","Dispute Resolution Only","Dispute Resolution Only"),
  estimate = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
  ci.lower = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
  ci.upper = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)  )   

estimates.s.na <- data.frame(
  test = c( "ACLED", "ACLED", "ACLED",
            "UCDP", "UCDP", "UCDP",
            "ACLED", "ACLED", 
            "UCDP", "UCDP"),
  model = c("One treatment", "Both treatments", "Incl. controls",
            "One treatment", "Both treatments", "Incl. controls",
            "Excl. controls", "Incl. controls",
            "Excl. controls", "Incl. controls"),
  variable = c("No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions",
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions",
               "No. of Subgroups \n w/ Policing Institutions", "No. of Subgroups \n w/ Policing Institutions",
               "No. of Subgroups x \n Constitutional Regulation","No. of Subgroups x \n Constitutional Regulation",
               "No. of Subgroups x \n Constitutional Regulation","No. of Subgroups x \n Constitutional Regulation" ),
  specification = c("Security Only", "Security Only","Security Only",
                    "Security Only", "Security Only","Security Only",
                    "Security Only", "Security Only","Security Only","Security Only"),
  estimate = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
  ci.lower = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
  ci.upper = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)  )   


estimates.robust <- bind_rows(
  list(
    estimates.robust,
    estimates.s.na,
    estimates.dr.na
  )
)


# Labels for plot
estimates.robust$specification[estimates.robust$specification == "Controlling for Jurisdictional Conflict"] <- "Controlling for \n Jurisdictional Conflict"
estimates.robust$specification[estimates.robust$specification == "Controlling for All Supgroups"] <- "Controlling for \n All Subgroups"
estimates.robust$variable[estimates.robust$variable == "No. of Subgroups x \n Constitutional Regulation (First Principal Component)"] <- "No. of Subgroups x \n Constitutional Regulation"
estimates.robust$variable[estimates.robust$variable == "Supergroup Policing x \n Constitutional Regulation (First Principal Component)"] <- "Supergroup Policing x \n Constitutional Regulation"
estimates.robust$variable[estimates.robust$variable == "No. of Subgroups x \n Constitutional Regulation (min-max)"] <- "No. of Subgroups x \n Constitutional Regulation"
estimates.robust$variable[estimates.robust$variable == "Supergroup Policing x \n Constitutional Regulation (min-max)"] <- "Supergroup Policing x \n Constitutional Regulation"

# Order
estimates.robust$specification_factor <- factor(estimates.robust$specification,
                                                levels = c(
                                                  "Constitutional Regulation: PC1",
                                                  "Constitutional Regulation: Min-Max",
                                                  "Non-Geo Outcome",
                                                  "Authority Conflicts Only",
                                                  "Territory Conflicts Only",
                                                  "Controlling for \n Jurisdictional Conflict",
                                                  "Controlling for \n All Subgroups",
                                                  "Security Only",
                                                  "Dispute Resolution Only",
                                                  "2015-2020",
                                                  "Groups outside \n Sub-Saharan Africa",
                                                  "All Groups",
                                                  "Mixed Effect", 
                                                  "Main"))

estimates.robust$model_factor <- factor(estimates.robust$model, 
                                        levels = c("Incl. controls", "Excl. controls", "Both treatments", "One treatment"))


# Plots group level
est_sup <- estimates.robust  %>%
  filter(variable=="Supergroup Policing") %>%
  ggplot(aes(x=estimate, y=specification_factor, shape=model_factor)) +
  geom_point(position = position_dodge(width=.3))+ 
  geom_errorbarh(aes(xmin=ci.lower, xmax=ci.upper), position=position_dodge(width=.3),
                 height=0) + 
  labs(x="Supergroup Policing Capacity", y="", shape="Specification:") + 
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_grid(. ~ test, scales = "free_x") +
  theme(panel.spacing = unit(2, "lines")) +
  scale_shape_discrete(breaks=c("One treatment", "Both treatments", "Incl. controls")) 


est_sub <- estimates.robust  %>%
  filter(variable=="No. of Subgroups \n w/ Policing Institutions") %>%
  ggplot(aes(x=estimate, y=specification_factor, shape=model_factor)) +
  geom_point(position = position_dodge(width=.3))+ 
  geom_errorbarh(aes(xmin=ci.lower, xmax=ci.upper), position=position_dodge(width=.3),
                 height=0) + 
  labs(x="No. of Subgroups w/ Policing Institutions", y="", shape="Specification:") + 
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_grid(. ~ test, scales = "free_x") +
  theme(panel.spacing = unit(2, "lines")) +
  scale_shape_discrete(breaks=c("One treatment", "Both treatments", "Incl. controls")) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "bottom") 

coef_robust_grp <- (est_sup + est_sub )  +
  plot_layout(guides = "collect") & theme(legend.position = 'bottom') 

ggsave("coef_robust_grp.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 25,
       height =16,
       units = "cm",
       dpi=700)


# Plots interaction
est_sup_reg <- estimates.robust  %>%
  filter(variable=="Supergroup Policing x \n Constitutional Regulation") %>%
  ggplot(aes(x=estimate, y=specification_factor, shape=model_factor)) +
  geom_point(position = position_dodge(width=.3))+ 
  geom_errorbarh(aes(xmin=ci.lower, xmax=ci.upper), position=position_dodge(width=.3),
                 height=0) + 
  labs(x="Supergroup Policing x Constitutional Regulation", y="", shape="Specification:") + 
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_grid(. ~ test, scales = "free_x") +
  theme(panel.spacing = unit(2, "lines")) +
  scale_shape_discrete(breaks=c("Excl. controls", "Incl. controls")) 

est_sub_reg <- estimates.robust  %>%
  filter(variable=="No. of Subgroups x \n Constitutional Regulation") %>%
  ggplot(aes(x=estimate, y=specification_factor, shape=model_factor)) +
  geom_point(position = position_dodge(width=.3))+ 
  geom_errorbarh(aes(xmin=ci.lower, xmax=ci.upper), position=position_dodge(width=.3),
                 height=0) + 
  labs(x="No. of Subgroups x Constitutional Regulation", y="", shape="Specification:") + 
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_grid(. ~ test, scales = "free_x") +
  theme(panel.spacing = unit(2, "lines")) +
  scale_shape_discrete(breaks=c("Excl. controls", "Incl. controls")) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "bottom")

coef_robust_reg <- (est_sup_reg + est_sub_reg )  +
  plot_layout(guides = "collect") & theme(legend.position = 'bottom') 

ggsave("coef_robust_reg.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 25,
       height =16,
       units = "cm",
       dpi=700)



##### Appendix I / Figure A11: Group-Level Interaction -----

#Note: functions specified above (Analysis Set-Up)

#Additional tratments and functions need to be specified for

  # Mixed effects
formulae.me.5 <- sapply(outcomes_lags, 
                        function(x) reformulate(c(treat.inter, x[2], baseline.me, controls.me), response = x[1]))
  # Dispute Resolution only
treat.inter.dr <- "disp_capacity*policing_subgroups"
formulae$twfe$twfe.dr.5 <- lapply(outcomes_lags, 
                                  construct_formula, treat = treat.inter.dr, baseline = baseline.twfe, controls = controls.twfe)
  #Security only
treat.inter.s <- "sec_capacity*policing_subgroups"
formulae$twfe$twfe.s.5 <- lapply(outcomes_lags, 
                                 construct_formula, treat = treat.inter.s, baseline = baseline.twfe, controls = controls.twfe)


five.ls <- list(
  lm.acled.5 <- feols(formulae$twfe$twfe.5[[1]],
                      data = nc_tpi, vcov = "hetero"), #estimates
  lm.ucdp.5 <- feols(formulae$twfe$twfe.5[[2]], 
                     data = nc_tpi, vcov = "hetero"), #estimates
  terr.lm.5 <- feols(formulae$twfe$twfe.5[[3]],
                     data = nc_tpi, vcov = "hetero"), #estimates.terr
  aut.lm.5 <- feols(formulae$twfe$twfe.5[[4]],
                    data = nc_tpi, vcov = "hetero"), #estimates.aut
  ng.lm.ucdp.5 <- feols(formulae$twfe$twfe.5[[9]],
                        data = nc_tpi, vcov = "hetero"), #estimates.ng
  future.lm.acled.5 <- feols(formulae$twfe$twfe.5[[1]],
                             data = nc_tpi[nc_tpi$year > 2014,]), #estimates.future
  future.lm.ucdp.5 <- feols(formulae$twfe$twfe.5[[2]],
                            data = nc_tpi[nc_tpi$year > 2014,]), #estimates.future
  all_grps.lm.acled.5 <- feols(formulae$twfe$twfe.5[[10]],
                               data = nc_rep, vcov = "hetero"), #estimates.all_grps
  all_grps.lm.ucdp.5 <- feols(formulae$twfe$twfe.5[[11]],
                              data = nc_rep, vcov = "hetero"), #estimates.all_grps
  nonAfr.lm.acled.5 <- feols(formulae$twfe$twfe.5[[1]],
                             data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]), #estimates.nonAfr
  nonAfr.lm.ucdp.5 <- feols(formulae$twfe$twfe.5[[2]],
                            data = nc_tpi[nc_tpi$region!="Sub-Saharan Africa",]), #estimates.nonAfr
  amar.lm.acled.5 <- feols(formulae$twfe$twfe.5[[7]],
                           data = nc_tpi, vcov = "hetero"), #estimates.subct
  amar.lm.ucdp.5 <- feols(formulae$twfe$twfe.5[[8]],
                          data = nc_tpi, vcov = "hetero"), #estimates.subct
  lmer.ucdp.5 <- lmer(formulae.me.5[[2]], 
                      data = nc_tpi), #estimates.lmer
  lmer.acled.5 <- lmer(formulae.me.5[[1]], 
                       data = nc_tpi), #estimates.lmer
  med.lm.ucdp.5 <- feols(formulae$twfe$twfe.5[[6]],
                         data = nc_tpi, vcov = "hetero"), #estimates.med
  med.lm.acled.5 <- feols(formulae$twfe$twfe.5[[5]],
                          data = nc_tpi, vcov = "hetero"), #estimates.med
  dr.lm.acled.5 <- feols(formulae$twfe$twfe.dr.5[[1]],
                         data = nc_tpi, vcov = "hetero"), #estimates.dr
  dr.lm.ucdp.5 <- feols(formulae$twfe$twfe.dr.5[[2]],
                        data = nc_tpi, vcov = "hetero"), #estimates.dr
  s.lm.acled.5 <- feols(formulae$twfe$twfe.s.5[[1]],
                        data = nc_tpi, vcov = "hetero"),  #estimates.s
  s.lm.ucdp.5 <- feols(formulae$twfe$twfe.s.5[[2]],
                       data = nc_tpi, vcov = "hetero") #estimates.s
)


# Extract estimates
estimates.five <- data.frame(
  test = c( "ACLED", "UCDP", "UCDP", "UCDP", "UCDP", 
            "ACLED", "UCDP", "ACLED", "UCDP", "ACLED", "UCDP",
            "ACLED", "UCDP", "ACLED", "UCDP", "ACLED", "UCDP", 
            "ACLED", "UCDP", "ACLED", "UCDP"),
  model = c("Incl. controls"),
  variable = c("Policing Capacity x \n Subgroup Count"),
  specification = c("Main", "Main", "Territory Conflicts Only", "Authority Conflicts Only",
                    "Non-Geo Outcome", "2015−2020", "2015−2020", "All Groups", "All Groups",
                    "Groups outside \n Sub-Saharan Africa", "Groups outside \n Sub-Saharan Africa",
                    "Controlling for \n All Subgroups", "Controlling for \n All Subgroups",
                    "Mixed Effect", "Mixed Effect",
                    "Controlling for \n Jurisdictional Conflict", "Controlling for \n Jurisdictional Conflict",
                    "Dispute Resolution Only", "Dispute Resolution Only",
                    "Security Only", "Security Only"),
  estimate = c(coef(lm.acled.5)["policing_capacity:policing_subgroups"],
               coef(lm.ucdp.5)["policing_capacity:policing_subgroups"],
               coef(terr.lm.5)["policing_capacity:policing_subgroups"],
               coef(aut.lm.5)["policing_capacity:policing_subgroups"],
               coef(ng.lm.ucdp.5)["policing_capacity:policing_subgroups"],
               coef(future.lm.acled.5)["policing_capacity:policing_subgroups"],
               coef(future.lm.ucdp.5)["policing_capacity:policing_subgroups"],
               coef(all_grps.lm.acled.5)["policing_capacity:policing_subgroups"],
               coef(all_grps.lm.ucdp.5)["policing_capacity:policing_subgroups"],
               coef(nonAfr.lm.acled.5)["policing_capacity:policing_subgroups"],
               coef(nonAfr.lm.ucdp.5)["policing_capacity:policing_subgroups"],
               coef(amar.lm.acled.5)["policing_capacity:policing_subgroups"],
               coef(amar.lm.ucdp.5)["policing_capacity:policing_subgroups"],
               coef(lmer.acled.5)$cow[1,"policing_capacity:policing_subgroups"],
               coef(lmer.ucdp.5)$cow[1,"policing_capacity:policing_subgroups"],
               coef(med.lm.acled.5)["policing_capacity:policing_subgroups"],
               coef(med.lm.ucdp.5)["policing_capacity:policing_subgroups"],
               coef(dr.lm.acled.5)["disp_capacity:policing_subgroups"],
               coef(dr.lm.ucdp.5)["disp_capacity:policing_subgroups"],
               coef(s.lm.acled.5)["sec_capacity:policing_subgroups"],
               coef(s.lm.ucdp.5)["sec_capacity:policing_subgroups"]),
  ci.lower = c(confint(lm.acled.5)["policing_capacity:policing_subgroups", 1],
               confint(lm.ucdp.5)["policing_capacity:policing_subgroups", 1],
               confint(terr.lm.5)["policing_capacity:policing_subgroups", 1],
               confint(aut.lm.5)["policing_capacity:policing_subgroups", 1],
               confint(ng.lm.ucdp.5)["policing_capacity:policing_subgroups", 1],
               confint(future.lm.acled.5)["policing_capacity:policing_subgroups", 1],
               confint(future.lm.ucdp.5)["policing_capacity:policing_subgroups", 1],
               confint(all_grps.lm.acled.5)["policing_capacity:policing_subgroups", 1],
               confint(all_grps.lm.ucdp.5)["policing_capacity:policing_subgroups", 1],
               confint(nonAfr.lm.acled.5)["policing_capacity:policing_subgroups", 1],
               confint(nonAfr.lm.ucdp.5)["policing_capacity:policing_subgroups", 1],
               confint(amar.lm.acled.5)["policing_capacity:policing_subgroups", 1],
               confint(amar.lm.ucdp.5)["policing_capacity:policing_subgroups", 1],
               confint(lmer.acled.5)["policing_capacity:policing_subgroups", 1],
               confint(lmer.ucdp.5)["policing_capacity:policing_subgroups", 1],
               confint(med.lm.acled.5)["policing_capacity:policing_subgroups", 1],
               confint(med.lm.ucdp.5)["policing_capacity:policing_subgroups", 1],
               confint(dr.lm.acled.5)["disp_capacity:policing_subgroups", 1],
               confint(dr.lm.ucdp.5)["disp_capacity:policing_subgroups", 1],
               confint(s.lm.acled.5)["sec_capacity:policing_subgroups", 1],
               confint(s.lm.ucdp.5)["sec_capacity:policing_subgroups", 1]),
  ci.upper = c(confint(lm.acled.5)["policing_capacity:policing_subgroups", 2],
               confint(lm.ucdp.5)["policing_capacity:policing_subgroups", 2],
               confint(terr.lm.5)["policing_capacity:policing_subgroups", 2],
               confint(aut.lm.5)["policing_capacity:policing_subgroups", 2],
               confint(ng.lm.ucdp.5)["policing_capacity:policing_subgroups", 2],
               confint(future.lm.acled.5)["policing_capacity:policing_subgroups", 2],
               confint(future.lm.ucdp.5)["policing_capacity:policing_subgroups", 2],
               confint(all_grps.lm.acled.5)["policing_capacity:policing_subgroups", 2],
               confint(all_grps.lm.ucdp.5)["policing_capacity:policing_subgroups", 2],
               confint(nonAfr.lm.acled.5)["policing_capacity:policing_subgroups", 2],
               confint(nonAfr.lm.ucdp.5)["policing_capacity:policing_subgroups", 2],
               confint(amar.lm.acled.5)["policing_capacity:policing_subgroups", 2],
               confint(amar.lm.ucdp.5)["policing_capacity:policing_subgroups", 2],
               confint(lmer.acled.5)["policing_capacity:policing_subgroups", 2],
               confint(lmer.ucdp.5)["policing_capacity:policing_subgroups", 2],
               confint(med.lm.acled.5)["policing_capacity:policing_subgroups", 2],
               confint(med.lm.ucdp.5)["policing_capacity:policing_subgroups", 2],
               confint(dr.lm.acled.5)["disp_capacity:policing_subgroups", 2],
               confint(dr.lm.ucdp.5)["disp_capacity:policing_subgroups", 2],
               confint(s.lm.acled.5)["sec_capacity:policing_subgroups", 2],
               confint(s.lm.ucdp.5)["sec_capacity:policing_subgroups", 2]))

save(estimates.five, file = "estimates/estimates.five.RData")

# Order
estimates.five$specification_factor <- factor(estimates.five$specification,
                                              levels = c(
                                                "Non-Geo Outcome",
                                                "Authority Conflicts Only",
                                                "Territory Conflicts Only",
                                                "Controlling for \n Jurisdictional Conflict",
                                                "Controlling for \n All Subgroups",
                                                "Security Only",
                                                "Dispute Resolution Only",
                                                "2015−2020",
                                                "Groups outside \n Sub-Saharan Africa",
                                                "All Groups",
                                                "Mixed Effect", 
                                                "Main"))


estimates.five %>%
  ggplot(aes(x=estimate, y=specification_factor)) +
  geom_point(position = position_dodge(width=.3))+ 
  geom_errorbarh(aes(xmin=ci.lower, xmax=ci.upper), position=position_dodge(width=.3),
                 height=0) + 
  labs(x="", y="") + 
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_grid(. ~ test, scales = "free_x") +
  theme(panel.spacing = unit(2, "lines"))

ggsave("coef_plot_grp_inter.pdf",
       plot = last_plot(),
       path = fig.path,
       device = "pdf",
       width = 18,
       height =10,
       units = "cm",
       dpi=700)





